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Abstract 

Least-squares Petrov-Galerkin (LSPG) model-reduction techniques such as the Gauss-Newton with Ap¬ 
proximated Tensors (GNAT) method have shown promise, as they have generated stable, accurate solutions 
for large-scale turbulent, compressible flow problems where standard Galerkin techniques have failed. How¬ 
ever, there has been limited comparative analysis of the two approaches. This is due in part to difficulties 
arising from the fact that Galerkin techniques perform optimal projection associated with residual minimiza¬ 
tion at the time-continuous level, while LSPG techniques do so at the time-discrete level. 

This work provides a detailed theoretical and computational comparison of the two techniques for two 
common classes of time integrators: linear multistep schemes and Runge-Kutta schemes. We present a 
number of new findings, including conditions under which the LSPG ROM has a time-continuous represen¬ 
tation, conditions under which the two techniques are equivalent, and time-discrete error bounds for the 
two approaches. Perhaps most surprisingly, we demonstrate both theoretically and computationally that 
decreasing the time step does not necessarily decrease the error for the LSPG ROM; instead, the time step 
should be ‘matched’ to the spectral content of the reduced basis. In numerical experiments carried out on a 
turbulent compressible-flow problem with over one million unknowns, we show that increasing the time step 
to an intermediate value decreases both the error and the simulation time of the LSPG reduced-order model 
by an order of magnitude. 

Keywords: model reduction, GNAT, least-squares Petrov-Galerkin projection, Galerkin projection, CFD 


1. Introduction 

While modeling and simulation of parameterized systems has become an essential tool in many indus¬ 
tries, the computational cost of executing high-fidelity simulations is infeasibly high for many time-critical 
applications. For example, real-time scenarios (e.g., model predictive control) require simulations to exe¬ 
cute in seconds or minutes, while many-query scenarios (e.g., statistical inversion) can require thousands of 
simulations corresponding to different parameter instances of the system. 

Reduced-order models (ROMs) have been developed to mitigate this computational bottleneck. First, 
they execute an offline stage during which computationally expensive training tasks (e.g., evaluating the 
high-fidelity model at several points in the parameter space) compute a representative low-dimensional 
‘trial’ basis for the system state. Then, during the inexpensive online stage, these methods quickly compute 
approximate solutions for arbitrary points in the parameter space via projection: they compute solutions 
in the span of the trial basis while enforcing the high-fidelity-model residual to be orthogonal to a low¬ 
dimensional ‘test’ basis. They also introduce other approximations in the presence of general nonlinearities 
or non-affine parameter dependence. See Ref. m and references within for a survey of current methods. 
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By far the most popular model-reduction technique for nonlinear ordinary differential equations (ODEs) 
is Galerkin projection [S3], wherein the test basis is set to be equal to the trial basis, which is often computed 
via proper orthogonal decomposition (POD) [35]. Galerkin projection can be considered continuous optimal , 
as the Galerkin-ROM velocity minimizes the ODE (time-continuous) residual in the £ 2 -norm. In addition, 
for specialized dynamical systems (e.g., Lagrangian dynamical systems), performing Galerkin projection is 
necessary to preserve problem structure [391 [3211331 . However, theoretical analyses—in the form of time- 
continuous error bounds [45] and stability analysis [3Tj —as well as numerical experiments have shown that 
Galerkin projection can lead to significant problems when applied to general nonlinear ODEs: instability 
[46] . inaccurate long-time responses [52] 03], and no guarantee of a priori convergence (i.e., adding basis 
vectors can degrade the solution) 071 Section 5]. In turbulent fluid flows, some of this poor performance can 
be attributed to the trial basis ‘filtering out’ small-scale modes essential for energy dissipation. 

To address these shortcomings, alternative projection techniques have been developed, particularly in 
fluid dynamics. These include stabilizing inner products mmm ; introducing dissipation via closure 
models il m M, nil ED] or numerical dissipation [ 36] : performing nonlinear Galerkin projection based on 
approximate inertial manifolds mmm-, including a pressure-term representation 03< [33] ; modifying 
the POD basis by including many modes (such that dissipative modes are captured), changing the norm 
[35], enabling adaptivity El, or including basis functions that resolve a range of scales [5] or respect the 
attractor’s power balance M: and performing Petrov-Galerkin projection [23] . 

Alternatively, a promising new model-reduction methodology eschews Galerkin projection in favor of 
performing projection at the fully discrete level , i.e., after the ODE has been discretized in time [191 . This 
discrete-optimal method -known as least-squares Petrov-Galerkin (LSPG) projection—computes the solu¬ 
tion that minimizes the £ 2 -norm of the (time-discrete) residual arising at each time step; this ensures that 
adding basis vectors yields a monotonic decrease in the least-squares objective function. When equipped 
with gappy POD [27] (a least-squares generalization of—and precursor to—the discrete empirical interpola¬ 
tion method [24] ) to approximate the discrete residual as a complexity-reduction mechanism, this approach 
is known as the Gauss-Newton with Approximated Tensors (GNAT) method [3T]. While LSPG projec¬ 
tion does not necessarily guarantee a priori accuracy and stablility for turbulent, compressible flows, it has 
been computationally demonstrated to generate accurate and stable responses for such problems on which 
Galerkin projection yielded unstable responses |19l f20l l2l|. 

In spite of its promise, theoretical analysis has been limited to developing consistency conditions for 
snapshot collection num and discrete-time error bounds for simple time integrators mm- In particular, 
major outstanding questions include: (1) What are time-continuous and time-discrete representations of the 
Galerkin and LSPG ROMs for broad classes of time integrators? (2) Are there conditions under which the 
two techniques are equivalent? (3) What discrete-time error bounds are available for the two techniques 
for broad classes time integrators? Related to the third issue is how parameters (e.g., time step or basis 
dimension) for the LSPG ROM should be chosen to optimize performance. This work aims to fill this gap 
by performing a number of detailed theoretical and computational studies that compare Galerkin and LSPG 
ROMs for the two most important classes of time integrators: linear multistep methods and Runge-Kutta 
schemes. We summarize the most important theoretical results (which map to the three questions above) as 
follows: 
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Figure 1: Relationship between Galerkin and discrete-optimal ROMs at the time-continuous (ODE) and time-discrete (OAE) 
levels. Bolded outlines imply the ROM associates with a minimum-residual solution at that time-discretization level. Dashed 
lines imply the relationships are valid under certain conditions (see Theorems |4.2| and |4.3| ). 
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1. Continuous and discrete representations 

• Galerkin projection and time discretization are commutative (Theorem |3.4[ ). 

• LSPG ROMs can be derived for Runge-Kutta schemes (Section |4~T] ). 

• The LSPG ROM has a time-continuous (i.e., ODE) representation under certain conditions (Sec¬ 
tion [472} Figure [l]). This ODE depends on the time step At. 

2. Equivalence conditions 


• Galerkin and LSPG ROMs are equivalent for explicit time integrators (Corollaries |5.1| and 5.21. 

• Galerkin and LSPG ROMs are equivalent in the limit of At —> 0 for implicit time integrators 
(Theorem |5.3[ ). 

• Galerkin ROMs are discrete optimal for symmetric-positive-definite residual Jacobians (Theorems 


5.4-5.61. 


3. Error analysis 


We provide local and global a posteriori and a priori error bounds for both Galerkin and LSPG 


ROMs for linear multistep schemes (Section 6.11. 

For backward differentiation formulas, we show that the LSPG ROM can yield a lower local a 
posteriori error bound than the Galerkin ROM (Corollary |6.4[ ) because it solves a time-global op¬ 
timization problem (over a time window) rather than a time-local optimization problem (Remark 


6.9). 


For the backward Euler time integrator, we show that an intermediate time step should yield the 


lowest error bound (Corollary 6.17 and Remark 6.18). 

• For the backward Euler time integrator, we show that a larger basis size leads to a smaller optimal 
time step for the LSPG ROM (Corollary 6.171. 

• We provide global a posteriori and a priori error bounds for both Galerkin and LSPG ROMs for 
Runge-Kutta schemes (Section [673] ) . 

Figure [l] summarizes time-continuous and time-discrete representations of the two techniques. 

In addition to the above theoretical results, we present numerical results for a large-scale compressible 
fluid-dynamics problem with turbulence modeling characterized by over one million degrees of freedom. 
These results illustrate the practical significance of the above theoretical results. Critically, we show that 
employing an intermediate time step for the LSPG ROM can decrease both the error and the simulation time 
by an order of magnitude, which is a highly non-intuitive result that is of immense practical significance. 

The remainder of the paper is organized as follows. Section [2] formulates the full-order model, including 
its representation at the time-continuous and time-discrete levels. Section [3] presents the Galerkin ROM 
at the continuous and discrete levels, and Section [4] does so for the LSPG ROM. In particular, we provide 
conditions under which the LSPG ROM has a time-continuous representation. Section [5] provides conditions 
under which the Galerkin and LSPG ROMs are equivalent; in particular, equivalence holds for explicit 


integrators (Section 5.1), in the limit of the time step going to zero for implicit integrators (Section 5.2), 
and for symmetric-positive-definite residual Jacobians (Section |5.3|). Section [6] provides error analysis for 
Galerkin and LSPG ROMs for linear multistep schemes (Section 6.1), Runge-Kutta schemes (Section 6.3), 
and a detailed analysis in the case of backward Euler (Section 6.2). Section [ 7 ] provides detailed numerical 
examples that illustrate the practical importance of the analysis. Finally, Section [8] provides conclusions. 

In the remainder of this paper, matrices are denoted by capitalized bold letters, vectors by lowercase bold 
letters, scalars by unbolded letters. The columns of a matrix A £ R mxn are denoted by Gq £ R m , i £ N(n) 
with N(a) := {1,..., a} such that A := [ai • • • a„]. We also define No(a) := {0,..., a}. The scalar-valued 
matrix elements are denoted by a.y £ M such that a,j := [a±j ■ ■ ■ a mj ] , j £ N(n). A superscript denotes 
the value of a variable at that time instance, e.g., x n is the value of x at time nAt, where At is the time 
step. 
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2. Full-order model 


We begin by formulating both the time-continuous (ODE) and time-discrete (OAE) representations of 
the full-order model (FOM). 

2.1. Continuous representation 

In this work, we consider the full-order model (FOM) to be an initial value problem characterized by a 
system of nonlinear ODEs 


x(0)=x 0: (2.1) 

where x : [0, T] —> denotes the (time-dependent) state, x 0 £ R^ denotes the initial condition, and 
/ : x [0,T] —> M. n with (£,r) i —> /(£, r). This ODE can arise, for example, from applying spatial 

discretization (e.g., finite element, finite volume, or finite difference) to a partial differential equation with 
time dependence. We note that most model-reduction techniques are applied to parameterized systems 
wherein the velocity f is also parameter dependent. However, we limit our presentation to unparameterized 
systems for notational simplicity, as we are interested comparing Galerkin and LSPG ROMs for a given 
instance of the ODE. 


2.2. Discrete representation 

A time-discretization method is required to solve Eq. (2.11 numerically. We now characterize the full- 
order-model OAE, which is the time-discrete representation of the model, for two classes of time integrators: 
linear multistep schemes and Runge-Kutta schemes. 


2.2.1. Linear multistep schemes 


A linear fc-step method applied to numerically solve Eq. (2.1) can be expressed as 

k 


j=0 


k 

aiX n ~i = At'^Pjf 
j=o 


,n-J j-n-j 


( 2 . 2 ) 


where At is the time step, the coefficients ay and /3j define a specific linear multistep scheme, ao ^ 0, and 
k 

22 a j = 0 is necessary for consistency. In this case, the OAE is characterized by the following system of 
i=o 

algebraic equations to be solved at each time instance n £ N(T/Af): 


r n (w n ) = 0, (2.3) 

where w n £ M. N is the unknown variable and r n : R^ —> denotes the linear multistep residual defined as 


r n ( w ) := aow — Atf3 0 f(w, t n ) + ^ ay 

i=i 

Then, the state can be updated explicitly as 


At £&/( 

3 = 1 


x n -Ct n - J 


(2.4) 


Hence, the unknown is equal to the state. These methods are implicit if /3 q 0. 
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2.2.2. Runge-Kutta schemes 

For an s-stage Runge-Kutta scheme, the OAE is characterized by the following system of algebraic 
equations to be solved at each time step: 


r? <) = 0, *eN(a). 

Here, the Runge-Kutta residual is defined as 

S 

rf (w 1: .. .,w s ) := Wi - /(a;" -1 + A t'^a ij w j ,t n ~ 1 + c z At ), i £ N(s) 

3 = 1 


and the state is explicitly updated as 

S 

x n =x n ~ 1 + At J2 b i w i- 

1=1 


(2.5) 


( 2 . 6 ) 


(2.7) 


Here, the unknowns correspond to the velocity dx/dt at times t n ~ l +c, At, i £ N(s), and the coefficients 
bi , Cj, and ay define a specific Runge-Kutta scheme. These methods are explicit if ay = 0, V) > i and are 
diagonally implicit if ay = 0, V) > i. 


3. Galerkin ROM 

This section provides the time-continuous and time-discrete representations of the Galerkin ROM, as well 
as key results related to optimality and commutativity of projection and time discretization. 


3.1. Continuous representation 

Galerkin-projection reduced-order models compute an approximate solution x sa x with x £ to 
Eq. ( |2.1[ ) by introducing two approximations. First, they restrict the approximate solution to lie in a low¬ 
dimensional affine trial subspace Xq + Ran (<&), where $ £ K JVxp with <I> = I denotes the given reduced 
basis (in matrix form) of dimension p -C TV and Ran (A) denotes the range of matrix A. This basis can 
be computed by a variety of techniques, e.g., eigenmode analysis, POD [22, or the reduced-basis method 
HUSHl [561S2U55]. Then, the approximate solution can be expressed as 

x(t) = x 0 + H>x(t), (3.1) 


where x : [0, T] —> R p denotes the (time-dependent) generalized coordinates of the approximate solution. 
Second, these methods substitute x ■£- x into Eq. (2.1) and enforce the ODE residual to be orthogonal to 
Ran(3>), which results in a low-dimensional system of nonlinear ODEs 

= $ T f( x o + t), ®(0) = 0. (3.2) 


Remark 3.1. In order to obtain computational efficiency, it is necessary to reduce the computational com¬ 
plexity of repeatedly computing matrix-vector products of the form 3? T /. This can be done using a variety 
of methods, e.g., collocation (303110], gappy POD [231101131101101], the discrete empirical interpolation 
method (DEIM) [12l [231221 EH [6] , reduced-order quadrature [5] , finite-element subassembly methods [3 [29] , 
or reduced-basis-sparsihcation techniques (23). However, in this work we limit ourselves to comparatively 
analyzing different projection techniques. For this reason, we do not perform additional analysis for such 
complexity-reduction mechanisms; this is the subject of follow-on work. 


We now restate the well-known result that Galerkin projection leads to a notion of minimum-residual 
optimality at the continuous level. This is reflected in the top-right box of Figure[l] where the bolded outline 
indicates minimum-residual optimality. 


= 


Theorem 3.2 (Galerkin: continuous optimality). If the reduced basis is orthogonal, i.e., 
then the Galerkin ROM (3.1)-(3.2) is continuous optimal in the sense that the approximated velocity mini¬ 
mizes theI 2 -norm of the FOM ODE residual (2.1) over Ran (<fr), i.e., 


dx 

— (x 0 + &x,t)=a,rg min 

dt i;6Ran(3») 


v - f(x 0 + 


(3.3) 
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Proof. Because % = 4>^jr, problem (3.3) can be written as 


da: 

dt 


(x 0 + $x,t) 


arg min g ( v ), 
*eR p 


(3.4) 


where g ( v ) := ||3>h — f(x 0 + 3>i:,f)|||. We now assess whether Eq. ( |3.4| ) holds, i.e., whether ^ as defined 
by Eq. ( |3.2| ) is the minimizer of g. 

The function g can be expressed as g ( v ) = v T ^v — 2v 1 & T f(x o + 3 > se , t)+f(xo+&x, t) T f(xo+<&x, t). 
Due to the strict convexity of the function g , the global minimizer v* is equal to the stationary point of g, 
i.e., v* satisfies 


0=^ (v*) = 2‘f> T $v* - 2<f> T f(x 0 + $x, t) 
dv v 7 

i,* = $ T f(x 0 + <f>x,t), 


(3.5) 

(3.6) 


where orthogonality of $ has been used. Comparing Eqs. (3.6) and (3.2) shows %( x 0 + &x,t) = v* 
is the desired result. □ 


which 


Remark 3.3 (Galerkin ROM enrichment yields a monotonic decrease in the FOM ODE residual). 

Due to optimality property (3.3) of the Galerkin ROM, adding vectors to the trial basis—which enriches 
the trial subspace Ran (3?)—results in a monotonic decrease in the minimum-residual objective function in 
problem (3.3), which is simply the £ 2 -norm of the FOM ODE residual. Because the FOM ODE residual is 
equivalent to the difference between the ROM and FOM velocities, this implies that the £ 2 -norm of the error- 
in the ROM velocity will nronotonically decrease as the trial subspace is enriched. 


Thus, Galerkin ROMs exhibit desirable properties (i.e., minimum-residual optimality) at the time- 
continuous level. We now derive a time-discrete representation for the Galerkin ROM, noting that these 
properties are lost at the time-discrete level. 


3.2. Discrete representation 


As before, a time-discretization method is needed to numerically solve Eq. (3.2). 
the OAE for the Galerkin ROM. 


We now characterize 


3.2.1. Linear multistep schemes 


A linear fc-step method applied to numerically solve Eq. (3.2) can be expressed as 

k k 

” / n-j +n-j 


djx 1 j = Atf 


3 = 0 


3=0 


Here, the OAE is characterized by the following system of algebraic equations to be solved at each time step: 

r n ( w n ) = 0. (3.7) 

Here, the discrete residual corresponds to 

k k 

r n ( w) := a 0 w — At(3o& T f(x 0 + t n ) + ^ a.jX n ~^ — At ^ f + i ^x n ~\f n ~ : '^j (3.8) 

3 =i i=i 

and the generalized state is explicitly updated as 
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3.2.2. Runge-Kutta schemes 


Applying an s-stage Runge-Kutta method to solve Eq. (|3.2|) leads to an OAE characterized by the 

(3.9) 


following system of algebraic equations to be solved at each time step: 

r n i =0, ieN(s). 
Here, discrete the residual is defined as 


(*i 


:= Wi — $> T f(x 0 + <&x n 1 + At ^ ciij&Wj , t n 1 +c i At), i G N(s) 


and the generalized state is computed explicitly as 


' 71 — 1 


3 = 1 


+ At bjw'l. 

i-l 


(3.10) 


(3.11) 


Note that the Galerkin-ROM solution satisfying Eqs. (3.7) or (3.9) does not in general associate with the 


solution to an optimization problem; therefore, the optimality property the method exhibits at the continuous 
level has been lost at the discrete level. We now show that Galerkin projection and time discretization are 
commutative; this implies that Galerkin ROMs can be analyzed, implemented, and interpreted equivalently 
at both the time-discrete and time-continuous levels. This corresponds to the rightmost part of Figure [l] 

Theorem 3.4 (Galerkin: commutativity of projection and time discretization). 

Performing a Galerkin projection on the governing ODE and subsequently applying time discretization yields 
the same model as first applying time discretization on the governing ODE and subsequently performing 
Galerkin projection. 


Proof. Linear multistep schemes. Eq. (3.7) was derived by performing Galerkin projection on the continuous 


representation of the FOM and subsequently applying time discretization. If instead we apply Galerkin 


projection to the discrete representation of the FOM in Eq. (2.3), set w = Xq + and x l = Xq + , 

oq = 0 and 4? T 4> = 7, we obtain the following OAE to be solved at each time step: 


i G N 

4? T r’ 


(n), and use J2 J= i' 

1 (xq + &w) = 0. Comparing the definition of the linear multistep residual (2.4) with Eq. (3.8) reveals 


r n (w) = 4? T r ra (xq + &w) , 


(3.12) 


which shows that the same discrete equations r n (w) = 0 are obtained at each time step regardless of the 
ordering of time discretization and Galerkin projection. 


Runge-Kutta schemes. Eq. (3.9) was derived by performing Galerkin projection on the continuous FOM 


representation and then applying time discretization. If instead we apply Galerkin projection to the discrete 
FOM representation in Eq. (2.5), set x n ~ x = Xq + Wi = i G N(s), and use 4> T 4> = I 


$ T r 


obtain the following OAE to be solved at each time step: 
the definition of the Runge-Kutta residual (|2.6[) with Eq. (|3.10| reveals 


Wi = <i'w„ * 

(4>thi,.... $w s ) = 0, 


i G N(5 


(tbi,..., w s ) = 4> T r" ($lt!i,..., $* s ), i G 


we 

Comparing 

(3.13) 


which shows that the same discrete equations r" (wi, ..., w s ) = 0, i G N(s) are obtained at each time step 
regardless of the ordering of time discretization and Galerkin projection. □ 


4. Least-squares Petrov Galerkin ROM 

Rather than performing minimum-residual optimal projection on the full-order model ODE (i.e., at the 
continuous level), this can be executed on the full-order model OAE (i.e., at the discrete level). Doing so 
enables discrete optimality , which contrasts with the continuous optimality exhibited by Galerkin projection. 
In particular, we consider optimal projections that minimize the discrete residual(s) (in some weighted £ 2 - 
norm) arising at each time instance. 

We note that other residual-minimizing approaches have been developed in the case of linear HU and non¬ 
linear j40j steady-state problems, and space-time solutions [25] ■ In addition, a recently developed approach 
p)( has suggested L 1 minimization of the residual arising at each time instance for hyperbolic problems. 
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4-1. Discrete representation 

We begin by developing the time-discrete representation for the LSPG ROM for both linear multistep 
schemes and Runge-Kutta schemes. The latter is a novel contribution, as previous work has derived discrete- 
optimal LSPG ROMs only for linear multistep schemes mm- Optimality of this approach corresponds to 
the bolded bottom-left box of Figure [lj 


4-1.1. Linear multistep schemes 

As before with Galerkin projection, discrete-optimal ROMs compute solutions using two approximations. 
First, they restrict the approximate solution to lie in the same low-dimensional affine trial subspace x £ 


a;o + Ran($) as Galerkin methods; thus, the approximate solution can be written according to Eq. (3.1l. In 


the case of linear multistep schemes, the unknown at time step n is simply the state, i.e., w 11 = x n , which 
implies that w n = Xq + <&in”. Second, the discrete-optimal ROM substitutes in” ■£- w n into the OAE (2.3) 


and solves a minimization problem to ensure the approximate solution is optimal in a minimum-residual 
sense at the discrete level: 

\A{z)r n {z)\\l (4.1) 


or equivalently 


w = arg 


w = arg mm 
zeRp 


mm 

:G£Co+Ran(<I>)" 


| A (x 0 + <&z) r n (x 0 + $>z) 


(4.2) 


Here, A £ H& zxN with z < N is a weighting matrix that enables the definition of a weighted (semi)norm. 
Examples of such reduced-order models include the LSPG method presented in Refs. mmm ( A = i), 
LSPG with collocation (A = P with P consisting of selected rows of the identity matrix) [40], and the 
related GNAT method [TO] HI] (A = (P3> r ) + P with 3> r a basis for the residual and the superscript + 
denoting the Moore-Penrose pseudoinverse). 


Note that necessary conditions for the solution to Eq. (4.2) corresponds to stationary of the objective 


function in Eq. (4.2), i.e., the solution satisfies 


\l/"(in") T r" (a:o + <I>in”) = 0, 


(4.3) 


where the entries of S'” £ M. Nxp are 


ipij(w) =a mi (x o + <l>w ) da,nl+ </> fcj r” (a: 0 + $w)+ 


dw k 


dr ? 


(4.4) 


a m i(x o + &w)a m i(x 0 + ^w)^A-(x 0 + $in)()> fcj , i £ N(TV), j £ N(p), 


where a repeated index implies summation. Because Eq. (4.3) corresponds to a Petrov-Galerkin projection 


with trial subspace Ran (S') and test subspace Ran (’S'), the discrete-optimal projection can be referred to 
as a least-squares Petrov-Galerkin projection mm- 


4-1.2. Runge-Kutta schemes 

LSPG ROMs for Runge-Kutta schemes also approximate the solution according to Eq. (3.1). However, 
because the unknown at time step n and stage i is the velocity at an intermediate time point, i.e., in” = 
x (f" _1 + CjAf) for i £ N(s), we have in" = S>:c (t" _1 + c^At) for this case. Then, these techniques substitute 
in" ■£- tb” into the OAE (|2.5h and solve the following minimum-residual problem: 


(in”,..., in”) = arg min ^ \\A t (*i, • ■ •, z a ) r” (zi, ..., z s ) \ 

Zi,...,;Zi)<ERan ) s ^ 
i=l 


(4.5) 


or equivalently 


(in”,..., in”) = arg min V || A t ($zi,..., $z s ) r” ($Zi,..., $z s ) \\%. 


(zi,.. .,z s )eR pXs z 

i=1 


(4.6) 









Here, A, g 


zxN 


, i g N(s) with z < N are weighting matrices. As before, the solution to Eq. (4.6) 


corresponds to a stationary point of the objective function, i.e., it satisfies 


i=i 


”) T r” ($*”,..., $*") = 0, i = 1 ,..., s, 


(4.7) 


where entries of the test bases tF”- g 

L J 


Jxp , *, j g N(s) are 






w s ) =[Ai\ uk (&w 1,..., $m s ) 


• • • , $* g ) 


9[t 


[A]«fc($™i, • • •, < &w s )[A i \ um {$w 1 ,..., $ih s ) 


g[g 

<9[«b 




(4.8) 


where [-] t j denotes entry (i,j) of the argument. This again leads to a least-squares Petrov-Galerkin inter¬ 
pretation for the discrete-optimal ROM. 

In the explicit or diagonally implicit cases, we can consider an alternative notion of discrete optimality. 
Explicit Runge-Kutta schemes are characterized by a^- = 0, Vj > i, while diagonally implicit Runge-Kutta 
(DIRK) schemes [2] are characterized by a,;, = 0, Vj > i. In these cases, solutions w ", i g N(s) can be 
computed sequentially, i.e., 

<?ror) = °> * = !»•■•,« 

with 

2-1 

qi(w) :=w- /( x n ~ x + A tagw + A ,t n ~ 1 + CjAt), i = l,...,s. (4.9) 

3 =1 

We can then formulate the following sequence of optimization problems to compute discrete minimum- 
residual approximations: 

w? = arg min \\Ai(z)q™ {z)\\l, i = (4.10) 

z6Ran($) 

or equivalently 

w? = arg mm \\A i ($z)q'?($z)\\l, i = l,...,s. (4.11) 

Here, the associated Petrov-Galerkin projection is 

¥?«) T qr?(**?) = 0, * = (4.12) 


with test-basis entries of 

m ik M = [A^i&w)^[q^w) + [A i ] UJ ($th)[A i ]„ m ($th)^%^^^ fc , (4.13) 

dw m owe 

Note that in the explicit case, the LSPG ROM generally requires solving a system of nonlinear equations 
at each Runge-Kutta stage. Because dq^/dw = J, the system of equations is linear if A, are constant 
matrices, and only an explicit solution update is required if A, = I and 4> T 4? = I. 


Remark 4.1 (LSPG ROM enrichment yields a monotonic decrease in the FOM OAE residual). 

Due to optimality property (4.1) of the LSPG ROM, adding vectors to the trial basis—which enriches the 
trial subspace Ran (4>)—results in a monotonic decrease in the minimum-residual objective function in prob¬ 
lem (4.11, which is simply the weighted £ 2 -norm of the FOM OAE residual associated with linear multistep 
schemes. This result also holds for LSPG ROMs applied to Runge-Kutta schemes, as the computed solutions 
satisfy alternative optimality properties in the implicit (4.5) and explicit/diagonally implicit (4.10) cases. 


We now see the critical tradeoff between the Galerkin and LSPG ROMs: Galerkin ROMs exhibit continu¬ 
ous (minimum-residual) optimality, while LSPG ROMs exhibit discrete optimality. Without further analysis, 
it is unclear which of these attributes is preferable. Numerical experiments (Section [?| and supporting error 
analysis (Section [6]) will highlight the benefits of discrete optimality over continuous optimality in practice. 
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4- 2. Continuous representation 

Because the LSPG ROM introduces approximations at the discrete level, it is unclear whether it can be 
interpreted at the continuous level. In fact, it has not previously been shown that a continuous representation 
of the LSPG ROM even exists. We now show that an ODE representation of the LSPG ROM does indeed 
exist for both linear multistep schemes and Runge-Kutta schemes under certain conditions; however, the 
ODE depends on the time step used to define the LSPG ROM. This associates with the top-left section of 
the relationship diagram in Figure |T| 

Theorem 4.2 (LSPG ROM continuous representation: linear multistep schemes). 

The LSPG ROM for linear multistep integrators is equivalent to applying a Petrov-Galerkin projection to 
the ODE with test basis (in matrix form) 

V(x,t) = A T A(a 0 I-Atf3 0 ^(x 0 + $x : t)^ & (4.14) 


and subsequently applying time integration with a linear multistep scheme with time step At if A is a constant 
matrix and (at least) one of the following conditions holds: 

1. /3j = 0, j > 1 (e.ga single-step method), 

2. the velocity f is linear in the state, or 

3. /3q = 0 (i.e., explicit schemes). 


Proof. Applying Petrov-Galerkin projection to the full-order model ODE (2.1) using a trial subspace Xq 
Ran(4>) and test subspace Ran(SP) yields the following ODE 


'&(x,t) T $-~ = '&(x,t) T f(x 0 + *(0)=0, (4.15) 

at 

which can be written in standard form as 

fiqrt / \ — 1 

— = 4>J 'Sf(x,t) T f(x 0 + i(0)=0. (4.16) 


Case 1 Applying a linear multistep time integrator with the stated assumption of /3 ? = 0, j > 1 to numerically 
solve Eq. (4.16) results in the following discrete equations to be solved at each time instance: 


a 0 y n - At/3 0 (*{y n , t n ) T *) ' ¥(»", t n ) T f{x Q + *y n ,t n ) + £ a 3 x n ~^ = 0. (4.17) 

3 = 1 


Pre-multiplying by 'Sf(y n ,t n ) 1 4> yields discrete equations r n (y n ) = 0 with residual 


r n (w) := ao'®'('U’, t n ) T &w — Atf3o'&(w 1 t n ) T f{X q + &w, t n ) -f t n ) T &x n 

j =i 


(4.18) 


Comparing Eqs. (4.18) and (2.4) reveals r n (w) = ^f(w,t n ) T r n (a:o + &w) and so the solution y n satisfies 

*(»", t n ) T r n (x 0 + &y n ) = 0. (4.19) 

Under the stated assumptions, we have dr n /dw(x) = a^I — At/3 0 ^(x,t n ) and so the LSPG test 
basis 41" defined in Eq. (4.41 is equal to the test basis in Eq. (4.14| evaluated at time instance n, i.e., 
^ n (w) = &(w,t n ). Therefore, the solution w n to the LSPG OAE (4.3) satisfies 

ty(w n , t n ) T r n (x 0 + <f?w n ) = 0. (4.20) 

This shows that w n = y n , i.e., the solutions to the LSPG OAE and the OAE obtained after applying 
Petrov-Galerkin projection with test basis SI f(x,t) defined by Eq. (4.14) to the full-order model ODE and 
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subsequently applying time integration are equivalent under the stated assumptions, which is the desired 
result. 

Case 2 In this case, the test basis is independent of the state, i.e., 


*(i) = A 1 A [ a 0 I - 


(4.21) 


Applying a linear multistep time integrator to solve Eq. (4.161 and subsequently pre-multiplying by the 
constant matrix 1 i'(t n ) T 3? yields the following discrete equations arising at each time step 


r n ( y n ) = 0, 

where the residual is defined as 

k 

:=a a 'S’(t n ) T ‘f>w - At/3 0 ^/{t n ) T f(x Q + $w,t n ) + ^ a j *(t n ) T &x n - j - 


(4.22) 


r (w) 


3 =1 


AC 


(4.23) 


1=i 


Comparing Eqs. (4.231 and (2.41 reveals r n (w) = 4 t(t n ) T r n (xq + Hrw) and so the solution y n satisfies 


V(t n ) T r n (xq + $y n ) = 0. 


(4.24) 


Under these assumptions, we have dr n /dw = aoI-At(3 0 df/d£(-, t n ) and so the LSPG test basis ’E" defined 
in Eq. (4.4) is equal to the test basis in Eq. ( 4.21| ) at time instance n, i.e., \f"(u;) = \H(t n ). Therefore, the 
LSPG OAE (4.3) can be expressed as 

’4 <(t n ) T r n (x 0 + <f>w n ) = 0. (4.25) 

This shows that w n = y n , i.e., the solutions to the LSPG OAE and the OAE obtained after applying 
Petrov-Galerkin projection with test basis 'U'(t) defined by Eq. (4.14) to the full-order model ODE and 
subsequently applying time integration are equivalent under the stated assumptions. 

Case 3 The assumption /3 0 = 0 results in a constant test basis 


= a 0 A A4>. 


(4.26) 


Applying a linear multistep time integrator to solve Eq. (4.16) and subsequently pre-multiplying by the 
constant matrix ’4 ,T 4> yields 

r n (: y n ) = 0, (4.27) 

which is to be solved at each time step with a residual defined as 

k k 

r n (w) := a 0 'Z' T &w - Atfo^fixo + ^w,t n ) + J^a^ T ^x n ~ j - AtJ^^A' T f (x 0 + ^x n ^ j ,t n ~ j ) ■ 


1=1 


1=1 


(4.28) 

As in Case 2, this leads to r" (w) = VlU r n (xq + $w). Because ^y(ai) = «o I, w e also again have ^ n (w) = 
S&. This leads to the desired result, as the OAEs for the LSPG ROM and the ROM obtained after applying 
Petrov-Galerkin projection with test basis '4' to the full-order model ODE and subsequently applying time 
integration both satisfy ’$? T r n (xo + $w n ) = 0 under the stated assumptions. □ 

We now provide conditions under which the LSPG ROM for Runge-Kutta schemes can be expressed as 
an ODE. 

Theorem 4.3 (LSPG ROM continuous representation: Runge—Kutta schemes). The LSPG ROM 
for Runge-Kutta integrators is equivalent to applying a Petrov-Galerkin projection to the ODE with test basis 
(in matrix form) 

— - aT A (t a - 














and subsequently applying time integration if Aj = A, V* are constant matrices and the integrator is a singly 
diagonally implicit Runge-Kutta (SDIRK) scheme, i.e., aij = 0, Vj > i and an = a, 

Proof. Applying Petrov-Galerkin projection to Eq. using a trial subspace x 0 + Ran (<&) and test 

subspace Ran(’i') yields the following ODE (in standard form) 


dx 

dt 


(&(x,t) T $^ '&{x,t) T f(x 0 + ®x,t), a(0) = 0. 


(4.30) 


Applying an SDIRK time integrator to numerically solve Eq. (4.30) results in the following sequence of 
discrete equations to be solved at each time step: 


V?~ [(^r 1 + At X] Oij y], r- 1 + c z Atf*}- 1 *^- 1 + AtJ2 OijVj, ^ + tkAtf 

3 =1 3 =1 

f(x 0 + <&x n 1 + At y ajj&yf, t n - 1 + c,At) = 0, i = 
j =i 

Pre-multiplying by + At J2 ]=i O'ijVjA 71-1 + CiAf) T <& yields the following discrete equations 

Qi (y?) =0; i = l,...,s 


(4.31) 


with residual 


i-1 

qi(w) :='I?(x n ~ 1 + Ataw + At'^ / aijy’j,t n ~ 1 + CjAf) T • 

3 =1 

i— 1 

&w — f{x o + + Ataw + At ^ aij&y™, t n ~ 1 + CjAf) 

i=i 


(4.32) 


Comparing Eqs. (4.32) and (4.9) reveals 


i -1 


<li(w) = ^{x n 1 + Atatb + At Y ajjy'fA 71 1 + CiAt) T qf ($>w) , i=l, ...,s 

j=i 


such that the solutions yf satisfy 


^(a: n 1 + At'^2a ij y],t n 1 + CiAt) T qf (yf) = 0, *=l,...,s. 

i=i 


Now, under the stated assumptions, we have 


i -1 


^p-(u) = I — Ata^j^(x n 1 + Atau + At^^aijWj ,t n 1 + c,At) 


j=i 


(4.33) 


such that the LSPG test basis 'E(' defined in Eq. (4.13) is related to the test basis in Eq. (4.29) as follows: 

i-1 

^^(w) = \I>( x n ~ 1 + Ataw + A t^2 a ijWj , t n ~ 1 + Cj At). 

j=i 


1 Note that summation on repeated indicies is not implied. 
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Therefore, the solutions to" to the LSPG OAE (4.12) satisfy 


'&(x r ' 


At ^2 ciijW ", 

3 = 1 


t n 1 + CiAt) T q" (ih”) = 0, i = l, ...,s. 


(4.34) 


Comparing Eqs. (4.33) and (4.34) shows that the to™ = y", i G N(s), i.e., the solutions to the LSPG 


OAE and the OAE obtained after applying Petrov-Galerkin projection with test basis 4 f(x,t) defined by 


Eq. (4.29) to the full-order model ODE and subsequently applying time integration are equivalent under the 


stated assumptions, which is the desired result. □ 

We now show that the LSPG ROM has a time-continuous representation for all explicit and single-state 
Runge-Kutta schemes, which include the widely used forward Euler, backward Euler, and implicit midpoint 
schemes. 


Corollary 4.4 (LSPG ROM continuous representation: explicit Runge Kutta). The LSPG ROM 
for Runge-Kutta integrators is equivalent to applying a Petrov-Galerkin projection to the ODE with test basis 
(in matrix form) 

4 f(x,t) = A t A$ 

and subsequently applying time integration if Ai = A, Vi are constant matrices and an explicit Runge-Kutta 
scheme is employed. 


Proof. Explicit Runge-Kutta schemes are characterized by a,j = 0, j > i and so they satisfy the conditions 
Theorem 14.31 with a = 0. □ 

Corollary 4.5 (LSPG ROM continuous representation: single-stage Runge-Kutta). The LSPG ROM 
for Runge-Kutta integrators is equivalent to applying a Petrov-Galerkin projection to the ODE with test basis 
(in matrix form) 

4>(x,f) = A t A ^1 - Atau^(x 0 + <&x,t)^J $ 

and subsequently applying time integration if Ai = A, Vi are constant matrices and a single-stage Runge- 
Kutta scheme is employed. 


Proof. Single-stage Runge-Kutta schemes are characterized by s = 1 and so they satisfy the conditions of 
Theorem |4.3| with a = an. □ 


5. Equivalence conditions 

This section performs theoretical analysis that highlights cases in which Galerkin and LSPG ROMs are 
equivalent, in which case the Galerkin ROM exhibits both continuous and discrete optimality. This suggests 
that time discretization and minimum-residual projection are commutative in these scenarios. Section |5.1| 
shows that equivalence holds for explicit time integrators, Section |5.2| demonstrates equivalence in the limit 
of At —► 0, and Section [5.3| shows equivalence in the case of symmetric-positive-definite residual Jacobians. 


5.1. Equivalence for explicit integrators 

Corollary 5.1 (Equivalence: explicit linear multistep scheme). 

Galerkin projection is equivalent to LSPG projection with A = J for explicit linear multistep schemes. 


Proof. In the case of e xplicit linear multistep schemes, 
Case 3 of Theorem 


4.2 


with A = —)=I, as 'S' = 4> in this case. 

%/ao ’ 


= 0 and so Galerkin projection corresponds to 

□ 


Corollary 5.2 (Equivalence: explicit Runge—Kutta scheme). Galerkin projection is equivalent to LSPG 
projection with A = I for explicit Runge-Kutta schemes. 


Proof. In the case of explicit Runge-Kutta schemes, a,,- = 0 V) > i and so Galerkin projection corresponds 
to Theorem |4.3| with A = I and a = 0, as 4/ = $ in this case. □ 
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5.2. Equivalence in the limit of At —► 0 

Theorem 5.3 (Limiting equivalence). 

In the limit of At —* 0, Galerkin projection is equivalent to LSPG projection with A = -A=/ for linear 
multistep schemes and A,, = I , i £ N(s) for Runge-Kutta schemes. 


Proof. Linear multistep schemes. Consider solving the LSPG OAE (4.3) with A = ^=J. Then, the test 
basis defined in Eq. (4.4) is simply 

1 dr n 

^ n (w) = --— (xq + &w) 

ao dw 


From Eq. (2.4), we can write the residual Jacobian as 

dr 


df, 


8w N =<*oI- At/3 0 g^(u,t n ). 


Therefore, we have 


1 


At—>-0 


At->0 Qo 


df, 


lim SP n (w) = lim — a 0 1 — At/3o-—(x 0 + &w,t n ) I $ = $ 


OH 


and so in the limit of At —> 0, the LSPG ROM solution satisfies 

lim (w) T r n (xq + <&w n ) = <& r r n (x$ + &w n ) = 0. 
At—>o v ' v ' 


(5.1) 


Because the Galerkin ROM solution also satisfies Eq. (5.1) (see Eq. (3.12) of Theorem |3.4| ), the two techniques 
are equivalent in this limit, which is the desired result. 


Runge-Kutta schemes. Consider solving the LSPG OAE (4.7) with A; = I, i £ N(s). Then, the test basis 
defined in Eq. (4.8) is simply 

(*!,•• -,Ws) = 


Now, from Eq. (|2.6|) the Jacobian can be expressed as 

df 

dH 


(tti,...,« s ) = ISij - Ataij^(x n 1 + A t'^2a i jU j ,t n 1 +c i At). 


dw 


3 = 1 


Therefore, we have 


lim = lim I ISij - Ata t j^(x n 1 + At y~] OjjUj, t n 1 + c*At) I $ = $6j. 


At—>-0 


At—>-0 


d£ 


3 =1 


and so in the limit of At —► 0, the LSPG ROM solution satisfies 


jim ^ ..., w s ) T r ” (<&•&",..., r” (x 0 + <&w n ) =0, i £ N(s). (5.2) 


j =i 


Because the Galerkin ROM solution also satisfies Eq. (5.2) (see Eq. (3.13) of Theorem 3.4), the two techniques 
are equivalent in this limit, which is the desired result. □ 
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5.3. Equivalence for symmetric-positive-definite residual Jacobians 

Theorem 5.4 (Equivalence: linear multistep schemes). In the case of linear multistep schemes, Galerkin 
projection is equivalent to LSPG projection with A(z) = U (z), where U is the Cholesky /octo^j of the 
residual-Jacobian inverse 

dr n 

dw 

if dr n /dw (w n ,t n ) = cxqI — Af/3 0 §| ( w n ,t n ) is symmetric positive definite and if 

= 0, Vz, k. (5.4) 

ow k 

Here, index notation has been used. 


= U T U, 


(5.3) 


Proof. Under the stated assumptions, the LSPG test basis defined in Eq. (4.41 is equal to the trial basis, 
i.e., 4t"(t(; ra ) = 4>. By invoking Eq. (3.12), we can see that the OAEs for the the LSPG ROM (4.3) and 
Galerkin ROM (|3.7| both satisfy < t > r n ( x 0 + = 0, which is the desired result. □ 


Theorem 5.5 (Equivalence: diagonally implicit Runge Kutta schemes). In the case of diagonally 
implicit Runge-Kutta schemes, Galerkin projection is equivalent to LSPG projection with Ai (z) = JJ_ i (z), 
where TJ_ i is the Cholesky factor of the residual-Jacobian inverse 

= UjUi , (5.5) 

if dq™/dw ( w n , t n ) = I — Atag^ ^ x n ~ 1 + A tauw n + At Y^j=\ aijiv™, t n_1 + CjAtj is symmetric positive 
definite and if 


dq " 

dw 


d[u t y 

div m 


tymk [q 


,nin 

i \ e 


= o, 


Vj, k. 


(5.6) 


Here, index notation has been used. 


Proof. Under the stated assumptions, the LSPG test basis defined in Eq. (|4.13| is equal to the trial basis, 
i.e., \&"(uj”) = $, i 6 N(s). By invoking Eq. (|3.13|), 

021) and Galerkin ROM Ol both satisfy $ T q? 


we can see that the OAEs for the the LSPG ROM 
= 0, i = 1 ,..., s, which is the desired result. □ 


Theorem 5.6 (Equivalence: implicit Runge—Kutta schemes). In the case of Runge-Kutta schemes, 
Galerkin projection exhibits discrete optimality if df n /dw ( w n ,t n ) is symmetric positive definite and if 


dug 

dwk 


<t>kjre = 0 , 


Vi, k. 


(5.7) 


Here, index notation has been used and U £ R sNxsN th e Cholesky factor of the residual-Jacobian inverse, 

i.e., 

= U T U. (5.8) 

Here, 


dr n 

dw 


Wi 


1 

CO 

s 

gr-H 

S* 

1_ 


$ 


£ R sN , r n :w^ 


e R sJV , 4? := 


W s 


1 

CO 

9 

S co 

s*. 

_1 


$ 


£ R sNxs P 


Its derivative can be computed by solving the Lyapunov equation -q^~ U + U 


dr n 1 1 d 2 r n dr n 1 ^ 

dw J dwdwfc |_ dw J 
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Proof. First, note that solution (m”,..., tb”) to the Galerkin OAE (3.13) equivalently satisfies 


<f> 7 r n ($■&") = 0, 


where w := 


. We are now precisely in the situation of Theorem 


solution is the solution to the (discrete) optimization problem 

minimize||t/ ( z)f n ( z ) 

2gRan($) 

under the assumed conditions. This objective function can be written equivalently as 


5.4 


the Galerkin 


(5.9) 


\\U ( z)r n (z)\\l = ^|| ^2 U ij ( z i’---’ z s)r](z 1 ,...,z s )\\Z, 

i=i j=i 

where tj i3 £ M. NxN denotes the ( i,j ) block of U. Therefore, the Galerkin ROM solution satisfies 


(5.10) 


arg min V H H Ai i • • • > ® z s) r] (&zi, •. ■, $z s ) \\%, 

f=1 j=1 


(5.11) 


where A,, = . Comparing Eqs. (5.11) and (4.6) reveals that the Galerkin ROM satisfies a slightly more 


general notion of discrete optimality than the LSPG schemes considered in this work. □ 

This analysis demonstrates that Galerkin projection exhibits discrete optimality when the residual Jaco¬ 
bian is symmetric positive definite. This is aligned with recent work that has shown Galerkin projection to 
be effective for Lagrangian dynamical systems [B2G21I23J—which are characterized by symmetric-positive- 
definite residual Jacobians—due to the fact that Galerkin projection preserves properties such as symplectic 
time evolution and energy conservation. For these reasons, using Galerkin projection is sensible for problems 
exhibiting these characteristics. 


6. Error analysis 

Ultimately, we are interested in assessing the state-space error between the (computed) time-discrete 
ROM solution and the (unknown) time-continuous FOM solution. This error comprises two contributions: 
the state-space error between (1) the time-continuous FOM and time-discrete FOM solutions (i.e., time- 
discretization error), and (2) the time-discrete ROM and time-discrete FOM solutions. This section focuses 
on the latter and performs time-discrete state-space error analyses for Galerkin and LSPG ROMs applied to 
different time integrators. 

Section [6~T| derives error bounds for the Galerkin and LSPG ROMs for linear multistep schemes. Here, 
Theorem |6.5| provides a posteriori bounds that depend on the ROM solution, while Theorem |6.11| reports 
a priori bounds. Section |6.2| provides a posteriori error bounds for the backward Euler scheme, as well as 
additional analyses that highlight the important role of the time step in the LSPG ROM, which is discussed 
in Remark |6.18| Section [6.3| derives ROM error bounds for Runge-Kutta schemes. Here, Theorem 6.19 
provides a posteriori bounds, Corollary |6.20| sp ecializes these results to explicit Runge-Kutta and DIRK 
schemes, and Theorem |6.21| and Corollary |6.23| report a priori bounds for Runge-Kutta schemes. 

6.1. Linear multistep schemes 

Here, we perform error analysis for implicit linear multistep schemes. We will use subscripts *, G and 
P to denote the solution to full-order model OAE (2.3), Galerkin ROM OAE (3.7), and the LSPG ROM 
OAE (4.3), respectively. We also acknowledge that linear multistep schemes with k > 1 usually employ 
different coefficients (3j and a 3 for different time instances; this is necessary because at time instance n, 
a maximum of n + 1 states is available from past history (starting with the initial condition at n = 0). 
Therefore, we allow for coefficients that depend on the time instance n, i.e., a™ and /3". In addition, we 
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define <P n := V n (xp) whose entries are defined by Eq. ( |4.4[ ). 
arising at each time instance n for linear multistep schemes as 


a o x * = PoAtf (xq + x™,t n ) + r" 


x 


n—k 




n— 1 


'n— k 




a n 0 x n G = fflAt* T f (x 0 + <f>x n G ,t n ) + r n G 
a n 0 x n P = ffiAt ((vI> n ) T *) (* n ) T f ( x 0 + 4>x n P , t n ) + f n P 


We can then write the discrete equations 




*2 = 0 

(6.1) 



x G = 0 

(6.2) 

~n-k 

Xp , • 

A 71 ” 1 
• > x P 

H> 

II 

o 

(6.3) 


where 


x n ~ k ,. 

.,*" _1 

:= f/3"A tf (x Q + x n - e ,t n ~ l ) - a^x n ~ e ) , 



l=i ' ' 


x n ~ k 7 . 

.,*- 1 

:= ]T (0?A t<S> T f * 0 + <f>x n -\t n - e ) - , 

(6.4) 


e= i ' ' 


x n ~ k ,. 


:= Y, (/3"At ((¥") r $) _1 (’®'") T / (*o + $x n ~ e ,t n -^ - a?x n -^ , 

i=\ ' ’ 



and x k := x k — Xq. We also dehne the FOM residuals at time instance n associated with the trajectories 
associated with the FOM, Galerkin ROM, and LSPG ROM OAEs as 


r?(x) ■- a^x - AtpQ f (*o + x) - f 
r G (x) := a%x - Atp%f(x 0 + x) - f 
r P (x) := a^x - Atffif(x 0 + x) - f” &x 
We dehne the Galerkin and LSPG operators as 


~,n-k n -1 

"*'* ) ■ ■ ■ ! ■*'* 

*xl-\...,*xl-' 


i—k 


<i>x, 


(6.5) 

( 6 . 6 ) 

(6.7) 


Y := $$ T , and P” := $ ((’4'”) T $) 1 (^") T , 
respectively, and Galerkin and LSPG state-space errors at time instance n as 

Sx G := — &x G , and dxfi := as" — 3>a;p, 

respectively. Note that the Galerkin operator V is an orthogonal projector (<J> is assumed to be orthogonal), 
while the LSPG operator P is an oblique projector. As the second argument in f does not play any role 
for linear multistep schemes (the time index always matches that of the first argument), will drop it for 
notational convenience in this section and in Section [6.21 


6.1.1. A posteriori error bounds 

We proceed by deriving a posteriori error bounds for the Galerkin and LSPG ROMs for linear multistep 
schemes. We assume Lipschitz continuity of f in the first argument: 

(Ai) There exists a constant n > 0 such that for x,y G 

|| f(x) - f{y)|| 2 < n\\x-y\\ 2 . 


Theorem 6.1 (Local a posteriori error bounds: linear multistep schemes). If (Ax) holds and At < 


at 


/( 


% «), Vj G 


(n), then 


k k 


ii^gII 2 <E £ " 

(I-Y)f (* 0 + *x n G - 1 ) 

+ E^ n 


0 

1=0 


2 £= 1 



k 


k 



l|S*p|| 2 <E e " 

(I - P”) f (* 0 + 

+E^ n 


0 

1=0 


2 i?=i 


z 

where we have defined ej 1 := A t/h m , 7 ™ := (a™ + \@T\ «A t)/h m , and h m := 

<|-|/3 0 ™| K At 


( 6 . 8 ) 

(6.9) 
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Proof. It is sufficient to show bound (6.261, as the arguments for (6.251 are similar. Let n be fixed but 
arbitrary, then subtracting Eq. (6.3) from Eq. (6.1) yields 


kiii<^pIi 2 < m At /(*o+ o - p"/ (*o+ **p) 


+ 


Sr™' 1 


( 6 . 10 ) 


where Srf, 1 := r™ [a 


,n-k n- 

'* > • ■ • ) 


- ®f r i 


-n—k *n— 

Cp , . . . , Xp 


. Adding and subtracting f (x 0 + &Xp) and 


applying the triangle inequality leads to 

K\\\Sx n P \\ 2 < \Pq\ At (J(7 - P") / (*o + *Xp) || 2 + ||/ (*o + <)-f (*o + **p) Q 
Invoking (Ai), and using At < \a^\ /|/3q | n, we deduce 

l^o I At 


Sr 


n— 1 


\\Sx™ P \\ 2 < 


(7 - P") f (x 0 + *x™ P ) 


1 


dr 


n—1 


Next, we will estimate 


Sr™- 1 


. Using the definition of r™, rf, from (6.4) we derive 


Sr 


n— 1 


< e m\ a< f (*o +*r‘) - p'7 (*o + 


+ 




. ( 6 . 11 ) 


( 6 . 12 ) 


(6.13) 


Adding and subtracting fix o + <l>Xp ( ), applying the triangle inequality in conjunction with (Ai) yields 


Sr™ p ~ 1 


K K 

< El^ n | At 7 - P ”) f { x ° + * x P~ e ) + E W «Ai +|a?|) 

<=i 2 &=i 


da; 


l — i 


(6.14) 


Then (6.12) and (6.14) imply Eq. (6.9) where we have used the definitions for e™ and 7 ". □ 

We now establish a result that aids interpretability, as it provides a connection between the terms in the 
a posteriori error bounds and the optimality properties of LSPG and Galerkin ROMs. First, we note that 


\(I -Y)f (x 0 + $XQ,t n ) || 2 = min \\y -f (x 0 + $x G , t n ) || 2 , 

y£Ran(cp) 


(6.15) 


due to the optimality property associated with orthogonal projection; this term appears in local a posteriori 
error bound ( 6 . 8 ). 

Lemma 6.2 (Oblique projection as discrete-residual minimization). If /3™ = 0, j > 1 (e.g., back¬ 
ward differentiation formulas), then 


\0S\ Ai||(I ^ v)/ (*o + *x n G ,t n ) || 2 = \\r™ G {*x n G )h 
\0S\ At||(/ - P ™)f (xo + *x™ P ,t n ) || 2 = \\r™p(*x n p)\\ 2 .. 

If additionally the LSPG ROM employs A = I, then 

\\ r p(*Xp)h = grin \\r n P {y)\\ 2 . 

y£ Ran(<P) 


Proof. From Eq. (6.3), we have 


05 AtP n f (®„ + *x n P , t ™) = E at<S>x n p - e - E Am> f (*o + $xp~ e , t n ~‘ 

e=o 1=1 


(6.16) 

(6.17) 

(6.18) 


(6.19) 


such that 

k k 

/3g At(7 - P™)f (x 0 + <f>xp, t™) = 0SAtf (® 0 + $x™ P , t™) - E a^x n P ~ e + Y P™AtPf (x 0 + $x n P ~ e , t™~^ 

( 6 . 20 ) 


e=o 


e=i 


K 

= —r n (x 0 + 3>x n p~ l ) - (7 - P) E 7”A tf (x 0 + t n ~ e ) 


( 6 . 21 ) 
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Then, if /3" = 0, I > 1, the final term vanishes and we have Eq. (6.171, where a similar result (i.e., Eq. (6.16)) 
holds for the Galerkin ROM. If the LSPG ROM employs A = I, then the LSPG ROM solution ip satisfies 
Eq. (4.2) with A = I, which yields Eq. (6.18). □ 

Corollary 6.3 (Discrete v. continuous residual minimization). If (if = 0, j > 1 (e.g., backward dif¬ 
ferentiation formulas), the LSPG ROM employs A = I, and xfP 1 = xff 1 , l = 1,..., k, then 


min ||r£(»)|| 2 = \ffi\ At\\(I - P")/ (*„ + t n ) || 2 

ySRan(#) 

< \0S\ At||(/ - V)/ (*o + t n ) || 2 = m At min \\y - / (x 0 + $x G , t n ) || 2 . 

j/GRan(*) 


( 6 . 22 ) 


Proof. Under the stated conditions, r G = rf and the optimality result Eq. (6.18) holds, yielding the desired 
result. □ 

Corollary 6.3 shows that discrete-residual minimization (i.e., LSPG projection) rather than continuous- 
residual minimization (i.e., Galerkin projection) can produce a smaller value of a term that appears in 
the a posteriori error bounds; for example, this term appears on the right-hand side of Eqs. (6.8)-(6.9). 
This can be interpreted as arising from the fact that discrete-residual minimization computes the LSPG 
ROM solution that minimizes the entire normed quantity, while continuous-residual minimization performs 
orthogonal projection of the velocity given the Galerkin ROM solution. 


Corollary 6.4 (LSPG can produce lower local a posteriori error bounds than Galerkin). Under 
the assumptions of Corollary \6.3\ and Theorem\6.1\ the local a posteriori error bound for the LSPG ROM 


(6.9) is smaller than that for the Galerkin ROM (6.8). 


Proof. Under the conditions of Theorem 6.1 Eqs. (6.9) and (6.8) are valid local a posteriori error bounds. 
If (3f = 0, j > 1, and xfT e = xfU 1 , i = 1,..., k, then these bounds simplify to 


\\Sx n G \\ 2 < 
\\6xf\\ 2 < 


\PS\At 


h m 

l/3o|A t 


(J-V)/(a 
(J - P n ) / ( 




<&X 


n—t 


k 

£=1 

k 

i—i 


Sx n P ~ l 


Sxfr 1 


(6.23) 

(6.24) 


respectively. Under the conditions of Corollary 6.3 inequality (6.22) holds and the right-hand side of (6.24) 
will be smaller than the right-hand side of (6.23). □ 

This is an interesting result, as it provides conditions under which the LSPG ROM produces a smaller 
local a posteriori error bound than the Galerkin ROM. It also provides some theoretical justification for 
the numerical experiments in Section [7] which use the three-point backward-differentiation formula, wherein 
the LSPG uniformly outperforms the Galerkin ROM. We now extend to results of Theorem 6.1 to obtain 
global a posteriori error bounds. 

Theorem 6.5 (Global a posteriori error bounds: linear multistep schemes). Under the assumptions 
of Theorem\6.1\ we have 


n—1 min (k,j) 

ii^gII 2 <E E 

j —o e=o 

n—1 min (k,j) 

ii^pii 2 <E E 

i=o e=o 


ifoiC? - 


!{0} (j - 


I fa) I 

E Ii< 

fa )eA(j-e) i—i 

lfa)l 

e rw 

(Vi)€A(j~t) i=l 


-Emil Ur, 


-Emil Vr, 


„n-j+e 


_n—j+t 


(7 - Y) / (aj 0 


$x n G -° 


J — p n-j+i 


(6.25) 


) f ( x o 


$x n - J 


2 

(6.26) 


Here, 1 a(x) denotes the indicator function, A(p) := {(r]i) \ rp £ N (k), = p}, and denotes the 

length of the tuple ( rji ). 
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Proof. Notice that the term 


i-: 


/ (®0 


<&Xp 0 


in inequality (6.9) corresponds to the error intro¬ 


duced at time instance i € N(n) from the state at time instance i — j with j € N(/c); this term always appears 
in the time-local error bound with coefficient £*■. Further, it contributes to the error at a given time instance 
n > i through appropriate products of 7 ™- For example, the product 7 i 72 _ 7i _3 provides one possible 
path for ‘traversing’ the time-local error bounds from time instance n to an earlier error contribution at 
time instance n — 4. Applying this notion more generally and using Sx° P = 0, the error can be bounded by 


induction according to inequality (6.26). □ 


The bounds in Theorem |6.5| can be considered a posteriori error bounds, as they depend on the ROM 
solutions xq and Xp and can thus be computed a posteriori if the Lipschitz constant k can be estimated. 
Note that the rightmost term in the Galerkin bound corresponds to the orthogonal projection error of / onto 
Ran(3>), while the LSPG bound entails an oblique projector that depends on the ROM solution. Because 
this oblique projection can associate with a discrete residual-minimization property (i.e., Corollary |6.3[ ), the 
LSPG ROM can yield smaller error bounds (i.e., Corollary 6.4). Also, the first term within square brackets 


corresponds to errors incurred at the current time instance n (i.e., via the leftmost term on the right-hand 
side of inequality (6.9)), while the second term corresponds to all possible ‘paths’ from current time instance 
n to the error contribution at previous time instances (i.e., the rightmost term on the right-hand side of 
inequality (6.9)). We also note the importance of the time-step condition At < a° 0 /(/3g k), as the stability 


constants appearing in bounds (6.251-(6.26) exhibit unbounded growth as the time step At approaches its 
upper limit, i.e., 

lim h k = 0, Vfc € N (n), 

At H“o|/(|/ 3 o | K ) 


which implies that 


lim 

A H a S|/(KI 


£g = 00 


and 


A M 


lim 


= 00, Vk € N(n). 


We now derive a simplified variant of these bounds. 


Theorem 6.6 (Simplified global a posteriori error bounds: linear multistep schemes). Under the 
assumptions of Theorem\6. 1\ if additionally At < |q:q |(1 — e) /(«|/3q I) 0 < e < 1. Then, 


||&EgII 2 < (k + l)|/3 max | At 


\\Sxp\\ 2 < (k + l)|/3 max | At 


Here, we have defined 



exp (t n ne 1 (\P*\/\a*\ + |/5q I/IQ’oI)) - 1 


max 


{k\a*\ - |aj$|) + {k\P*\ + |/3ol)«Af ieN („ 

'k\a*\ \ H ex P (* nKe_1 (l/ 3 *l/l«*l + \Po\/K\)) - 1 


max 


{k\a*\ - |ckq|) + (fc \j3*\ + |/3q|)«A t je N(n) 

teNo(^) 


l«ol -|/3oK Ai : = 

j£N(n) 


|a*| +1/3*1 nAt := max 

jeN(n), teN(fe) 


kA t 


at 


ft 


kA t 


\Pn 


max |/3»| 

j'eN(n), teNo(fc) 


k > (.* := max fp, 

jGN(n) 


:= argrnax £/ 

teNo(fe) 


(I - P") / (® 0 + <S>x n p ~ e ) 


(7 - Y ) / (® 0 + * 4 ) 

(6.27) 

^7 - P ^ f + <f>x 

(6.28) 

(6.29) 

(6.30) 

(6.31) 

(6.32) 


p 
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Proof. We proceed by proving bound (6.28); the proof for bound (6.27) is similar. First, we define 


£ n := argmax 7 " 

^eNo(fc) 


5x n P ~ t 


as well as a ‘path’ from time instance n backward to the initial time by defining £( 0 ) = n with 

£(j + l)=£(j)-e CU) , j = 0,..., n — 1, 


with n < n and C(n) = 0. Then, from local a posteriori error bound (6.9), we have 


\\5x n P \\ 2 <(k + l)el 


(I - P n ) f (x 0 + <f>x n p *") 


+ k^ 


5x n p~* 


(6.33) 


(6.34) 


(6.35) 


where := |/3 max | At/flogl — |/3 qI K ^t) and := (|a*| + |/3*| KAtj/flajjl — |/3 q| kA t). Via recursion, we then 
obtain 


n—l / j—1 

ii^n 2 < e kj n vt) 

j—0 \ m—0 

n—l 

< (fc + l)|/3 max |E^ 

i =0 


(k + 1 ) e fio) 




3>a; 




(6.36) 


At 


| a* | + |/3*|«At \ 

Kl - |/3 q|«A< J |aj$| - \/3$\KAt 


( j _ P £ 0)) / ^ 0 + 


<t>:r 




2 

(6.37) 


n—l 


< (fc + 1)|/3 n 


= (k + l)| A. 


< (fc + l)|/3 n 


V* z.m ( ) 

^ l|a5|-|^l«At) 


At 


, m —0 


max 

(/ - F') / (x 0 + Qx^) 

J j'GN(n) 

\ J \ J 


At 


At 


1 ,ni |a* 1 + 1/3* |«At \ n i 


k(\a*\ 
( k\a* 


max 


|/3*|acA t) - |a#| + |^o|KAt jeN(n) 

teN 0 (£*) 


^eN 0 (t*) 


(/ - F') / (. 


Xq + <f>X j p 1 


(6.38) 

(6.39) 


exp 


( t n ne- 1 (|/3* 


+ 1 A 


0 1/ l u 0 


d) 


- 1 


max 


{k\a*\ - |ckqI) + {k |/3*| + \/3$\)nAt je N(n) 

teNo(^*)' 


(l - / (x 0 + 

(6.40) 




where we have used 


|a*| + |/3*|/v.At\ _ / |a*| 
l«ol - l^oI kA V \\ a o\ 


1 + nAt\P*\/\o*\ 
1 - «At|/3^|/|a5| 


laS 


1 + 


KAt(\/3*\/\a*\ + IAqI/IQqI) 
1 - «At|/3g |/|aj| 


the relation (l+x) n < exp (nx), and the following result (with x = l+«At|/3*|/|a*| and y = 1 — acA t|/3 q |/1«q|); 
if x > y, then (x — y)/y < e~ 1 (x — y) if and only if y > e > 0 . □ 

We note that due to the optimality property associated with orthogonal projection, we can write bound 
(16.271) equivalently as 


\\6x n G \\ 2 <\f) B 


At 


'k\a*\\ n ex P ( rKe 1 (IA*I/KI + IAol/l a ol)) 


- 1 




,,, , , lx ^ n 4 max mm 

(fc|a*| - |q!q|) + (k\/3*\ + |/3q|)kA t jen(n) y e Ran(*) 


y - f (x 0 + $x 


di.ll) 

We now prove conditions under which the a posteriori error bound is independent of the time step At and 
total number of time instances n; the bound is fixed for a given time t n . 

Corollary 6.7 (Time-step-independent global a posteriori error bounds: linear multistep schemes). 

Under the assumptions of Theorem 6.6 if additionally k\a*\ = |aj| (e.g., backward Euler, where k = 1 and 
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|cc*| = |ckq| = 1 ) then the global a posteriori error bounds (6.27)-(6.28) are independent of the time step and 
simplify to 


||I r n|| / (fc + l)|/3 max | 
i|,5:EgI12 ^ 


ir nil < (fc + l)IAnax| 

P " 2 - (fc|/3*| + |/3 0 *|) K 


exp (t n Ke 1 (|/3*|/|a*| + |/3ol/|a$|)) - l) max (I-V) f (x 0 

\ / / jGN( n) \ 


exp 


(t n /ce _1 (|/3*|/|a*| + |/3 q l/l a ol)) 


— 1 ) max 

j'eN(n) 
^ 6 N 0 {£*) 


<t>x J c 


(6.42) 

(/ - F') / (x 0 + $x j P ~ e ) 

(6.43) 


Proof. The result can be derived by substituting k\a*\ = |a(j| into inequalities (6.27) and (6.28). □ 

We now prove conditions under which a posteriori error bounds (6.27) (6.28) can be written in ‘residual 
form’, i.e., in terms of the discrete residual arising at each time step. This will enable the respective optimality 
properties of the Galerkin and LSPG ROMs to be compared in Remark |6.9| 


Corollary 6.8 (Simplified global a posteriori error bounds in residual form: linear multistep schemes). 

Under the assumptions of Theorem 6.6 if additionally fdj 1 = 0, j > 1, Vm £ N(n) (e.g., backward differen¬ 
tiation formulas). Then, 


II**gII 2 <(* + i) 


I<5*pII 2 — (^ +1) 


k\c 


k\c 


‘exp (Vbce- 1 (|/3*|/|a*| + |/3 0 *|/KI)) 

(fc|a*| - |ctq|) + (k\P*\ + |/SoD^At 

exp (t n Ke 1 (|/3*|/|a*| + |So I/1I)) 


— 1 


max Wr^^x- 1 


f6N(n) 


Gdl2 


- 1 


{k\a*\ - |a$|) + (k |/3*| + \P^\)nAt je N(n 


max ||rp($* J 


p) 112- 


(6.44) 

(6.45) 


If additionally A = I, then 

'fc|a*|Y l ex p( inKe_1 G^* 


||5 ®£|| 2 <(* + !) 


IAj 


5D) 


- 1 


max mm 


07, 


(fc|o*| - |a(5|) + (k |/3*| + |^q|)kA t ieN(n) 3 /GRan(#) 


\ r p{y)\[ 


(6.46) 


Proof. Under the stated assumptions /3 max = (do, £% = 0, and Lemma 6.2 holds, yielding the desired result. 

□ 


Remark 6.9 (Optimality in a posteriori error bounds). Comparing inequalities (6.41) and (6.46) high¬ 
lights the differences in how optimality affects the Galerkin and LSPG ROM error bounds. Writing these 
expressions more compactly under the conditions of Corollary \6.8\ yields 


1(5 x 


Gll 2 ^ v/3 0 At max min \\y - f (x 0 + <&x J G 

jGN(n) i/GRan(<[>) V 


||<fe #|| 2 < v max min 


where 


v := (k + 1 ) 


ji£N(n) 2 6a;o+Ran(<l») 

' fc | a *|\ n exp (W " 1 (|/3*|/ 


r J p (z) 


(6.47) 

(6.48) 


m/K\)) 


- 1 


(fc|a*| - |o#|) + (k |/3*| + \^\)nAt 


6.1.2. A priori error bounds 

We now derive a priori error bounds by slightly modifying the steps in the above proofs. The most 
significant difference in the subsequent results is that the oblique projection associated with the LSPG 
ROM no longer associates with residual minimization, as the argument of the operators corresponds to the 
full-order-model solution. 
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£=0 


where we have defined Gfi := |/3™| A t/h m 7 ™ := (|a™| + |/3™| KAi||P m || 2 )//i m > h m := |a™| — l/^o^I ftAf||P m || 2 - 
Other quantities are defined in Theorem\6.5[ 


Proof. Instead of adding and subtracting f (xq + &x P ) between Eqs. (6.10) and (6.11) and between Eqs. (6.131 


and (6.14) in the proof of Theorem 6.5 adding and subtracting P”/ (xq + x™) and using ||V ||2 = 1 yields 


the stated result. □ 


Note that bound (6.50) is not quite an a priori bound, as the operator P” depends on the LSPG ROM 


solution xfi; while this dependence could be removed, the bound in its current form facilitates comparison 
with the Galerkin bound. The same dependence persists for the remaining a priori error bounds in this 
section. 

Corollary 6.11 (Global a priori error bounds: linear multistep schemes). Under the assumptions 
of Corollary 6.10\ we have 


7 i—l min(fej) 

IGdl . t 


i«ii 2 <E E 

!{o e n^; m=i??m 

n-j+l 

0 

11 

'Sp 

0 

II 

»= 1 


7 i—l min(kj’) 

IGdl 


ii^ii 2 <E E 

i{o}(j-^)+ e m=i??m 

-n-j+i 

0 

11 

'Sp 

0 

II 

(vi)eA(j-i) »=1 



(J — Y) / (, 


x 0 + x+ 


(6.51) 


(j _ P »-*+*) f (x 0 + *?" J ') 


2 

(6.52) 


Proof. The result can be derived by following the same steps as Theorem 6.5 based on the local bounds in 
Corollary |6.10| □ 

Corollary 6.12 (Simplified global a priori error bounds: linear multistep schemes). Under the as¬ 
sumptions of Corollary 6.1C\ if additionally At < |oq|(1 — c)/(/«|/?q|) with 0 < e < 1 for the Galerkin ROM , 
and At < |dol(l — e)/( k|/3q11|||Pq|| 2 ) with 0 < e < 1 for the LSPG ROM, then 


H&EGII 2 + 1 )|/3max| At 


'k\a*\\ U ex P 1 0^*l/l a *l + l/WKI)) - 1 


max 


(k\a*\ - |oq|) + (k\/3*\ + \/3£\)nAt jeN(n 


(I-Y)f(x o + xl) 
(6.53) 


||<5a3jp|| 2 <(k + l)|/3max| At 


ex P (*"« e 1 (\P*\\\ v *h/\ a *\ + I 11111 ^°o 11 2 /1 ctq I) ^ - 1 


(fc|d*| - \a* 0 \) + (fc|/3*|||P*|| 2 + |^|||PS|| 2 )«At 


max 

leN(n) 

^eNo(tJ) 


(/ - F*) / (*0 


(6.54) 


Corollary 6.10 (Local a priori error bounds: linear multistep schemes). //( Ai) holds, At < a 3 0 /(/3g n), 
Mj £ N(n) for the Galerkin ROM, and At < a 3 0 /(/3q K||P n ||), V) £ N(n) for the LSPG ROM, then 

k k 

II^gII 2 <E £ ” (I-y)f{xv+<~ £ ) +Y.^\ 5 x g~% (6-49) 

e=o 2 e -1 

k k 

||^|| 2 <E^ (I-V n )f(x 0 + x:- e ) + E7,1^-1L (6-50) 
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Here, we have defined 


l«ol | A} | « A ^ll F oll 2 : = mm 

jeN(n) 


|a*| +1/3*1 nAt\\¥*\\ 2 := max 

. jeN(n), ^eN(fe) 


KAt||F || 2 


iAt\\¥ j 


k > it ■= max £ J i, 
je N(n) 


it := argmax sf 
^eNo(fe) 


(j - p") f (x 0 +*r £ ) 


Proof. The result can be derived by following the same steps as Theorem 6.6 based on the local bounds in 
Corollary |6.10 □ 

We now demonstrate conditions under which the a priori error bound is independent of the time step At. 

Corollary 6.13 (Time-step-independent global a priori error bounds: linear multistep schemes). 

Under the assumptions of Corollary \6.12\ if additionally fc|a*| = |oq| (e.g., backward Euler) for the Galerkin 
ROM, and fc|cC| = jdJI (e.g., backward Euler) for the LSPG ROM, then the a priori error bounds (6.53)- 
(6.54) are independent of the time step and simplify to 


\\Sx^\\ 2 < 


\Sx n P \\ 2 < 


(k m + i/3 0 *i)k 

(fc + 1) |/?max| 

(fc|^|||p *|| 2 + |^|||Psi| 2 )« 


— 1 ) max 

jeN(ra) 


(k + 1)l/?maxl (exp (t n Ke _1 (|r|/|«*| + l/WKD) 

(exp (r Ke -- 1 (r|||P*|| 2 /|d*| + |A 


(/ — ¥)/ (x 0 + x{ 


0II2/I a 0 


- 1 


max 

j£N(n) 

^eNo(t*) 


(/ - F') / (®o + xt^ 


(6.55) 


(6.56) 


Proof. As in Corollary 6.7 the result can be shown by substitution of the appropriate assumptions (i.e., 
k\a*\ = |aj5| for the Galerkin ROM, /c|cC| = |d(5| for the LSPG ROM) into inequalities (6.53) and (6.54). □ 


Remark 6.14 (Optimality in a priori error bounds). Because the argument of f in the a priori error 
bounds corresponds to the full-order-model solution, it is not possible to relate such quantities to ROM 


residuals as was done in the case of a posteriori error bounds in Lemma 6.2 As such, it is not possible 
to associate the oblique projection of the LSPG ROM with minimizing any component of the a priori error 
bounds, as was shown for a posteriori error bounds in Corollary \ 6. ■/[ However, it is possible to associate 
Galerkin projection with minimizing terms in the a priori error bounds, i.e., the following optimality result 
holds under the conditions of Corollary\6.12\ 


ll^rlU < v max min 

j'eN(n) 3 /SRan($) 


y 


f (®o 


(6.57) 


This is analogous to inequality (6.47) for a posteriori error bounds, where 

exp ^PV.e -1 (|/3" 


v ■.= (k + l)|/3„ 



«* + 


l/WKO) 


- 1 


(k\a*\ - |o#|) + (k\(3*\ + |/3*|)acA t 


6.2. Backward Euler 

In practice, error bounds for a specific time integrator can be derived by substituting the appropriate 
values of a™ and into the appropriate error bounds.; doing so can lend additional insight into the error 
bound. We now perform this exercise for the backward Euler scheme—which is a single-step method that 
can be characterized by Eq. \2.2\ with k = 1, «o = 1, oq = —1, /3o = 1, and /3i = 0—and selected bounds. 

Thus, that result 


We first note that the backward Euler scheme satisfies /?? = 0, j > 1 in Corollary 6.4 


provides conditions under which the LSPG ROM has a lower local a posteriori error bound than the Galerkin 
ROM. Next, we specialize the global a posteriori error bounds in Theorem |6.5| to the case of the backward 
Euler scheme. 
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Corollary 6.15 ( a posteriori error bounds: backward Euler). Under the assumptions of Theorem 6.1 
for the backward Euler scheme we obtain 


n— 1 


ll^clla < At J2 JJNJ+i ( 7 - V ) / (*o + J ) 


3=0 
n— 1 


1 5x 


P II 2 


< At J2 


3=0 


(hy +i 


(hy +i 


(i - / (*o + *x n P - J ) 


(6.58) 

(6.59) 


where h := 1 — nAt. Note that the time-step condition corresponds to At < 1/n in this case. 


Proof. Because it is a single-step method, the linear multistep coefficients do not vary between time instances; 
as a result, the constants appearing in Theorem |6.5| are also time-step independent and are h = 1 — nAt, 


£o = A t/h, £i = 0, 70 = (nAt + l)/h, and 71 = 1/h. Substituting these values into error bound (6.261 and 
noting that 


A(j) = 


{( 1 ) •••>!)€ 


j > 0 

otherwise, 


yields 


\\Sx n p \\ 2 [!{0 }ti) + (Wl N (n-l)(j)] (i P"- J ) / ( 


n— 1 


3 = 0 '- 


x 0 + &x n p j 


(6.60) 


(i/h.y 


which simplifies to bound (6.59). Derivation of bound (6.58) is identical and is thus omitted. □ 

We now specialize the results of the time-step-independent global a priori error bounds in Corollary |6. 13 
to the backward Euler scheme. 

Corollary 6.16 (Time-step-independent global a priori error bounds: backward Euler). Under the 
assumptions of Corollary \6. lf\ — which can be satisfied by backward Euler as k\a*\ = |ckq| for this scheme—we 
obtain the following for the backward Euler scheme: 


nii ^ o exp (t n Ke *) - 1 
\ox%Wo < 2 - max 

K j£N(n) 


|(Sa:p|| 2 < 2 


exp(t 


n .~ 11 


■ 0 ll 2 


)-i 


(i-v)/(* 0 +®i) 

(/ - F') / (x 0 + ®i) 


max 

j£N(n] 


(6.61) 
(6.62) 

'6lla« .' ' 

Note that the time step conditions correspond to At < (1 — e)/n for the Galerkin ROM and At < (1 — 
e)/(«||IP q || 2 ) for the LSPG ROM in this case. 

Proof. For backward Euler, k = |a*| = |qtq| = |/3q| = 1 and £* = |/3*| = 0. Substituting these quantities into 
bounds (6.55) and (6.56) produces the desired result. □ 

We now derive a result that highlights the (surprising) role the time step At plays in the LSPG error 
bound. As will be shown in the numerical experiments in Section [TJ this theoretical results can have an 
important effect on the performance of the LSPG ROM in practice. 

Corollary 6.17. If x solves an auxiliary problem that computes the full-space solution increment centered 
on the LSPG ROM trajectory 


Atf (a;o + x j ^j + $x j p 1 , j e N(n), 


then the following holds: 


n— 1 


||^|| 2 <(l + «At)^^ T 


3=0 

n—1 _ r 


= Af(l + nAt) 


3=0 


( h) j+ 1 


ll/(^)l| 2 . 


(6.63) 

(6.64) 

(6.65) 
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Here, p? := 


<J>A Xp — Ax 3 


* i ~ 7 —1 

p *aj p 

h 3 '■= h J / 


and Ax 3 := x 3 

As- 7 1| 2 . 


denotes the difference in solution increments at time instance j, where Ax 3 p := 
— ‘fcccjr 1 . We denote the relative solution increment at time instance j by 


Proof. Eq. ( |6.59| ) in conjunction with (6.17) (with|/?Q | = 1 and the appropriate form of the discrete residual 
for backward Euler) implies 


n— 1 


ii^pII 2 < E 


3=0 


(hV +1 


$Aajp 3 


At/( 


x 0 


$ x n p ~ 3 


( 6 . 66 ) 


We can also write the auxiliary equation (6.63) as Ax 3 = A tf (ad). j £ N(n), which allows us to rewrite 
bound (16.661) as 


iifepii 2 < E 


3=0 


(h ) 3+1 


($A x n P ~ j - Ax n ~ 3 ) - 


At ( / ^a; 0 + 3>A x p 3 + <l>x 


n—j—1 

P 


- f (x 0 + Ax 


* n ~ 3 + 


Lipschitz continuity of / leads to the bound (6.64). To obtain Eq. (6.651, we multiply and divide by 
||Aic w- ^||2 for each term in the summation and use Ax n ~i = A tf (x n ~ 3 ). □ 

Corollary 6.17 is useful in that it expresses the LSPG ROM error in terms of the (time-local) single-step 
errors incurred by projection along the LSPG ROM trajectory. In addition, this result highlights the critical 
role of the time step At in the performance of the LSPG ROM; the following remark provides this discussion. 


Remark 6.18. The time step At in the error bound (6.65) for the LSPG ROM solution plays an important 
role. In particular, decreasing the time step produces both beneficial effects (bound decrease) and deleterious 
effects (bound increase), which we denote by and respectively as follows: 


+ The time-discretization error decreases (this does not appear in the time-discrete error analysis above). 
- The number of overall time instances n increases, so there are more terms in the summation. 

+ The terms At(l + nAt) and l/(hy +1 decrease. 

? The term p, n ~ J may increase or decrease, depending on the spectral content of the basis 3?. 


We now discuss this final ambiguous effect. The term fi n can be interpreted as the relative error in solution 
increment over [(n — 1)A£, nAt], Clearly, the ability of the LSPG ROM to make p, n small depends on the 
spectral content of the basis <!?: if the basis only captures modes that evolve over long time scales, then fi n 
will be large (i.e., close to one), as the basis does not contain the ‘fast evolving’ solution components that 
change over a single time step. This suggests that the time step should be ‘matched’ to the spectral content 
of the reduced basis 3?. In Section [73] of the experiments, we explore this issue numerically, and demonstrate 
that the error bound is minimized for an intermediate value of the time step Af. 

We note that the above arguments do not hold for the Galerkin ROM, which is simply an ODE that does 
not depend on the time step. Instead, decreasing the time step should increase accuracy, as it has the effect 
of reducing the time-discretization error. 


6.3. Runge-Kutta schemes 

We now derive Runge-Kutta error bounds for the Galerkin ROM (3.10) and the LSPG ROM (4.7). For 
notational simplicity, we define /" : £ /(£,t" -1 +CiAt), i £ N(s), n £ N(T/Af). Rather than present the 

full collection of results as was done for linear multistep schemes in Section |6.1[ for compactness we instead 
focus only on the most important results for Runge-Kutta schemes: global a posteriori and a priori error 
bounds, as well as time-step-independent a priori error bounds. 
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We rewrite Eqs. (2.5), (3.9), and (4.7) as 


S 

w +,i = fi (xo + x™- 1 + At E a ij w *,j) > ieN(s) 

i=i 

(6.67) 

S 

*G,i = $T /” (*o + $*G _1 + At E ) * e N(s) 

1=1 

(6.68) 

*P,i = ((^”) T ^)“ 1 (^”) T /”(*o + + At i) 

l=i 


-(^E*)” 1 E (*r e ) T ( $ *p, e - /”(*o + + At e a «i**p,i) ) 

e=l,e^2 ' j=l ^ 

, i G N(s) (6.69) 

and set x® = Xq = xf, = 0. 



6.3.1. A posteriori error bounds 

We first derive a posteriori error bounds. 

Theorem 6.19 (Global a posteriori error bounds: Runge Kutta schemes). If (Ax) holds and At 
is such that 

(a) the matrix D G R sxs with entries dij := Sij — nAt\a,ij\ is invertible, and 

(b) for every x, y > 0, if Dx < y then x < D~ 1 y, 
then 


n— 1 


£=0 


||<5®gI| 2 < + l ] ki 


k—1 i —1 


s s s 

(E N- V)/r*(*o + + A t^Oij^w. 


(6.70) 


k=l i—1 


n—t 

G,j 


1=1 


and 


\\8x$h < (l + «Ai E NE^ 

1 =0 k —1 i=l 

f E nE^e^c/ - o++ a*E“«**^) 

\fc=l i=l 1=1 

-EnI>-E *((*ry*) _1 E (®rr(^-/r'(*o+*4?r'- 1 +AtE^^ 


fe=l 2—1 


e=l,e^2 


j=l 


2 / 

(6.71) 


where P? := $ ((’4'”) T $) 1 (¥? \) T . Here, inequalities applied to vectors hold entry wise. 


Proof. Galerkin ROM . First we will show bound (6.701. Subtracting Eq. (6.68) from Eq. (6.671 and applying 
the triangle inequality yields 


l «, il |2 < ||/?(*0 


_1 + At E' 

3=1 


■Vf?(x 0 + &X 


n— 1 
G 


A'E' J "‘ ,> " ? f ,; ,/ 

3=1 


i G N(s), 
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where Swq := to" ,• — Adding and subtracting /" ^a; 0 + 1 + At J2j =1 a ij^'u } G,jJ and invoking 

assumption (Ai), we deduce 


Id w 


G,i I 


5=1 


,j|||<5u; 


<3,51 


< 


(J-vj/r^o+saS^+A* 


5=1 


'o' 1, "h';,, 


+ac|| <5cc' 


n— 1 1 
G 


i e N(s). 


Selecting At small enough such that (a) and ( b ) hold yields 


s s s 

||<5^, fc || 2 < || 2 +Ac||A a3 -- 1 || 2 ^[Z^- 1 ] fci , (6.72) 


2=1 


3 = 1 


2—1 


where [*]^ denotes entry (i,j) of the argument. From explicit state updates (2.7) and (3.11), we obtain 

s 

\\SxM 2 < \\Sx%-% + AtJ 2 Mll^lh- 


k=1 


Using upper bound (6.72) yields 


s s 

||<5^|| 2 <(l + K At^|6 fe |^[r>- 1 ] fei )||«5 ;i: ”- 1 || 2 

A:—1 2—1 

s s s 

+At e \b k \ y][£) _i ]fci|| (i - v)/r(* 0 +^r 1 +a tE^**^ 


k —1 i —1 


5=1 


Finally, an induction argument produces the desired result (6.70). 

LSPG ROM . We now prove bound (6.711. Subtracting (6.691 from (6.67) and applying the triangle inequality 
yields 


\\Sw n P A \ 2 < 


S 

fi(x 0 + < _1 + At^« 


5=i 


S 

??/?(*o + *4p _1 + A 


7 p,5 


5=1 


+ 


*((*S) T *) _1 ± + 


e=l,e^2 


3 = 1 


i € N(s). 


Adding and subtracting /™ ^aj 0 + ‘F&p 1 + At J]j=i a ij^'u } p,j) and invoking assumption (Ai), for i G N(s) 
we deduce 


S S 

||^|| 2 - K At^|^-|||(5^|| 2 < |(I-P”)/r(*o + ^r 1 + At^a y -$t 


7 P,5 


5=1 

+ 

e=l,e^2 

Again selecting At small enough such that (a) and (6) hold yields 


5=i 


■ k\\Sx 


P II 2 


*((*s ) T *) _1 E (n) T (^p,e-/:(*o + ^r 1 +AtE^^ 


5=1 


||<5u5 


P,fel|2 < 


s s s 

E[0 _1 ]w|(/- P?)/?(*0 + ^P _1 + AtE^-^Po) | + «II^P _1 ||2 E^ _1 ]« 


5=1 


+ E[£> _1 W *((*”) T *) 1 E (^ef (**P,e - /e (*0 + $&P -1 + At E 


e=l,e^2 


J=1 


(6.73) 
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The explicit state updates from (2.7) and (3.11) and the bound (6.73) yield 


3 =1 


dij&Wpj 


\\Sx n P \\ 2 <(l + |6 fc | E[£> _1 ]fci) II^p _1 ||2 

k—1 i—1 

s s 

+ At £ |6 fc | - p D/r(*o + 

fc —1 2 = 1 

s s 

+ Af^|6 fe |^[£>- 1 ] fez . 

fc = l 2=1 

$((^) T ^) _1 E (K) T (^,.-/?(*o+**r 1 +AtE:^^ 


e=l,e^i 


i=i 


An induction argument yields the bound ( 6.71| ). □ 

As with linear multistep schemes, the error bound for the Galerkin ROM in (6.70) depends on the 
orthogonal projection error of /” onto Ran (3?), while the LSPG ROM error bound depends on an oblique 
projection; however, because the oblique projector depends on the LSPG ROM solution, this bound can 
be smaller than the Galerkin bound (as was demonstrated in Corollary 6.4—note that backward Euler is 
also a Runge-Kutta scheme). Further, notice that the complex ‘path traversing’ that appears in the linear 
multistep error bounds is not present for Runge-Kutta schemes; this is a consequence of the fact that previous 
time instances only influence the error at the current time step through induction for Runge-Kutta schemes. 

In addition, note that the LSPG ROM error bound is more complex than the Galerkin bound; the final line 
in bound (6.71) is a consequence of the test-basis coupling (4.8) for general implicit Runge-Kutta schemes. 
Finally, we point out that both bounds grow exponentially in time due to the amplification factor. We now 
present a simpler version of this error bound for explicit Runge-Kutta and DIRK schemes. 

Corollary 6.20 (Global a posteriori LSPG ROM error bound: explicit Runge Kutta and DIRK). 

Under the assumptions of Theorem 6.19 for explicit RK (6 = 1) and DIRK (6 = 0) schemes, we have 


n— 1 


£=0 


||&Ep || 2 < A(l + kA^E \ b k\ E[-° 

k—1 i- 


\ki 


, k=l 


2=1 


[i-'. 


'~ l )/r* (*0 + Sip-' -1 + A tJ2 


n—£ 

J P,3 


3= 1 


(6.74) 


Proof. For explicit and diagonally implicit Runge-Kutta schemes (4.12), we have = 0 when i ^ j. The 
proof is then an immediate consequence of (6.71). □ 


Note that the LSPG error bound (6.74) resembles the Galerkin error bound (6.70) much more closely 
than the previous LSPG bound (6.71), as explicit Runge-Kutta and DIRK schemes remove the coupling of 
the test basis across stages. 

6.3.2. A priori error bounds 


We now state the a priori versions of the Galerkin Runge-Kutta schemes (6.70) and the LSPG Runge- 
Kutta schemes (6.71). 

Theorem 6.21 (Global a priori error bounds: Runge-Kutta schemes). If (Ai) holds and At is 
such that 


(a) the matrices D £ R sxs (defined in Theorem 6.19) and D 
/vAt|oj i7 -|||P "||2 are invertible, and 


l sxs with entries [D”].j 


Sij 


(b) for every x,y > 0, if Dx < y then x < D l y and if D m x < y then x < [D m ] 1 y, m £ N(n) , 
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then 


e=o 


\\6xM2 < (l + KAt^lhl El D_1 ] 

) k= 1 2=1 

S S 

(Ew | ( J - W (*o+*r <_1 + At e || 2 ) ■ 


fc=1 i=l 


1=1 


and 


ii^i| 2 <AtE n (i+^ei^iE^ 

£—0 m—O k—1 2=1 

f E i 6 * i E[[£ n ~V 1 ]«|| (i - pt*)/?-* (*o+^r^- 1 +At e a ijW :;/) 

\k=1 i=1 

+ Ei& fc iEP n 'V 1 


1=1 


/C2- 


fc = l 2=1 


*((*s~v*) _ E (n~r ^-/r^-o+^-^+AtEoe,^-- 

e=l,e^2 ^ j=l 


(6.75) 


2 / 

(6.76) 


where we have used the convention that the empty product is equal to one. 


Proof. The proof of (6.75) follows the Galerkin-ROM derivation in Theorem 6.19 wherein the quantity 


Y f™ (x 0 +xf 1 +AtJ2 S j=i a ij w *j) is added and subtracted rather than f™(xo+&XQ 1 +AtJf l S j=i a : 

ows the LSPG 
instead of /” ^a; 0 + &x 




,*o 


x. 


AtEU 


anto': 


and |j"V|| 2 = 1 is used. On the other hand, ( |6.76[ ) follows the LSPG-ROM derivation by adding and subtract¬ 
ing P ?/?(a 




AtEj=i a b$™Pj)- D 


Similarly to Corollary 6.16 we now derive a time-step-independent variant of the error bound (6.751. 


A similar result can be shown for the LSPG bound (6.76); however, we omit this result for simplicity and 


instead will provide (more readily interpretable) a priori LSPG ROM error bounds for explicit Runge-Kutta 
and DIRK schemes in Corollary |6.24| 

Corollary 6.22 (Time-step-independent global a priori Galerkin ROM error bound: RK schemes). 

Under the assumptions of Theorem 6.21 assume additionally At < (1 —ui)/(Ka*) with a* := 

0 < u> < 1. Then, 


Hj Ulj 


1 2 and 


where b* := max ie ^(s) 1^1- 


max 
-i) 

2GN(s) 


s 

(i - Y)rr e (*o+ ^r *- 1 +At e || 2 ) > 


1 = 1 


(6.77) 


Proof. As in Corollary 6.16 by using an upper bound for the right hand side in (6.75) we obtain 


n— i s s p s s 

wsxzh < At E ( l + E n E^ -1 ^) ■ (E n E D_1 


ki 


£=0 


k—1 2=1 


k=l 2=1 


II (J - v>fr‘(*o+<-’- 1 +At £ y 

i£N(s) 


1=i 


(6.78) 
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We now derive a bound for the term ^2 \bk\ ^2 [D as 


k =1 


E i&fciE[^ _1 ]fci < ^EE^Ew ^ &*EE \ i d ^u < < &** V2 ii d 


k =1 


/C —1 2=1 


fc=l 2=1 


< 


&*s 3 / 2 


b* s 3 / 2 


(6.79) 


0min(-/") 0"max(^At[|Gqjj]jj) 1 ftAtu* 


where we have used the vector-norm equivalence relation ||ai||i < \/7V||£cH 2 with x £ R N , the matrix-norm 
equivalence relation ||A|| F < \/]V||A ||2 with A £ R NxN , the relation ||^4. —1 1| 2 = l/<r m i n (A), er min (A + 
B) > c m i n (A) — cr max (B) = cr m i„(A) — |jS || 2 (with A = I and B = — |]ij))s an d the assumption 

At < (1 — w)/(ko*). 

We can now compute an upper bound on the constant in inequality (6.78|) as 


n— 1 s s g s s 

At e (1 + «At e Hi E^ -1 ]*) ■ (E nE^" 1 ^ 


e=o 


k =1 2=1 

n —1 


k=l 


A x ' / _ /AAA 

^ A *E 1 +WT 


e=o 


KAtb*s 3/2 


1 — nAta* 


i=l 


/ b*s 3 / 2 \ 

VI — kA ta* ) 


+ kA t(b*s 3 / 2 - a *)^ (1 - nAta*)~ n - 1 


(6.80) 


(6.81) 


Then, setting t n = nAt we obtain 


^1 + kA t(b*s 3/2 - a*)^ (1 - kA ta*)~ n < exp(f n Kw~ Vs 3/2 ) 


(6.82) 


As in Theorem 


6.1 


we have used the relation ( 1 +a;)" < exp(n:r), and the result (with x = 1 +KAtl b*s 3 ^ 2 —a > 


and y = 1 — kA ta*) that states if x > y, then (x — y)/y < to 1 (x — y) if and only if y > to > 0. Finally, 
substituting in (6.82) in (6.81) and combining the resulting expression with (6.78) yields the desired result. 
□ 

As with linear multistep schemes, the rightmost term in the Galerkin a priori bound will always be 
smaller than that for the LSPG bound, as the former associates with an orthogonal projection error of a 
fixed vector. In addition, the LSPG bound depends on the LSPG ROM solution; while this dependence could 
be removed, the bound in its current form facilitates comparison with the Galerkin bound. We again notice 
the complex structure of the estimator in (6.76) compared to (6.75). To better understand the behavior the 
LSPG estimator in (6.76) we consider two subcases: explicit Runge-Kutta and DIRK schemes. 

Corollary 6.23 (Global a priori LSPG ROM error bounds: explicit Runge Kutta and DIRK). 

Under the assumptions of Theorem 6.21 for explicit RK (6 = 1) and DIRK (6 = 0) schemes, we have 

71-1 i-i 


\\Sx n p \\ 2 < AtE n (l + KAfElMEl^ 

£—0 m —0 


l-li 


J ki 


k =1 2=1 

E n EK £,n "V]w| - p rVE(* 0+zr*- 1 +a* e <*<?) 

, k— 1 2=1 j =1 


(6.83) 


Proof. For explicit and diagonally implicit Runge-Kutta schemes (4.12), we have = 0 when i 7 ^ j. The 
proof is then an immediate consequence of (6.76). □ 

Owing to the fact that the explicit Runge-Kutta and DIRK schemes removes the coupling of the basis, 
we again notice that the LSPG error bound (6.83) resembles the Galerkin ROM error bound (6.75). 


31 


























Corollary 6.24 (Time-step-independent global a priori LSPG ROM error bound: explicit RK and DIRK). 

Under the assumptions of Theorem 6.21 for explicit RK ( 9 = 1) and DIRK (9 = 0) schemes, assume addi¬ 
tionally At < (1 — u>)/(na*) with a* := min£ e N 0 (n_i) \\M\nhu\ 2 and 0 < u) < 1. Then, 


Wteph < 


( exp(t n KU)~ 1 b*s 3 / 2 ) - 1 

V K 


ll (J - + At £ ID- 


max 
-i) 

6N(s) 


3=1 


(6.84) 


Proof. The proof follows closely that of Corollary |6.22| We begin by deriving an upper bound for the right 
hand side in (6.751 as 


71-1 £-1 


\\6x n P \\ 2 <AtY, n 1 + kA^IM^D 


n—m -1 _i ] 


■ r^n—Z-T — ii 


£—0 m —0 


k—1 i— 1 


E 

,k= 1 1=1 

(«££„ l< J - o++ a* e ««»;,-') y • 

ieN(s) 


(6.85) 


i =i 


Analogously to the derivation in Eq. (6.79), we bound the term Ylk=l IM £2i= i[[-D" ^ Hfci as 

6 *s 3 / 2 

[V \ “ 1 w < T 

fc=l i=l 


EnE^ 


1 — kA ta* 


( 6 . 86 ) 


Finally, we compute an upper bound on the constant in inequality (6.851 as 

TL —1 £—1 s s / s s \ 

A *E ■ Ew DP" - '] -1 ]*) < a^e (f + 

£=0 m—0 ... - ' ' 7 * ' 


n—1 


k—1 i—1 


k k—1 i—1 


£=0 


kA tb*s 3/2 \ / b*s 3/2 


1 — kA ta* 


/ V's"'* \ 

\ 1 — kA ta* ) 


As this result is identical to upper bound (6.801 in the proof of Corollary 6.22 with a* replacing a*, we obtain 
the desired result by applying the remaining steps as Corollary |6.22| □ 


7. Numerical experiments 


This section compares the performance of Galerkin and LSPG ROMs on a computational-fluid-dynamics 
(CFD) application using a basis constructed by proper orthogonal decomposition. These experiments high¬ 
light the importance of the previous analyses, in particular the limiting equivalence of Galerkin and LSPG 
ROMs (Theorem 5.3), superior accuracy of the LSPG ROM compared with the Galerkin ROM (Corol¬ 
lary 16.41) , and performance improvement of the LSPG ROM when an intermediate time step is selected 


(Corollary |6. 17 and Remark 6.18). 

Note that these experiments could be carried out on any dynamical system yielding a system of nonlinear 
ODEs ( |2.1| ); we have selected compressible turbulent fluid dynamics due to both its wide interest and 
challenging nature: limited progress has been made to date on developing robust, accurate ROMs for such 
problems. The numerical experiments highlight this fact, as standard Galerkin ROMs generate unstable 
responses in all cases. 


7.1. Problem description 

The Galerkin and LSPG ROMs are implemented in AERO-F [Ml 130) . a massively parallel compressible- 
flow solver. AERO-F solves the steady or unsteady compressible Navier-Stokes equations with various 
closure models available for turbulent flow, and employs a second-order node-centered finite-volume scheme. 
For model-reduction algorithms, all linear least-squares problems and singular value decompositions are 
computed in parallel using the ScaLAPACK library HU- 
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The full-order model corresponds to an unsteady Navier-Stokes simulation of a two-dimensional open 
cavity with a length-to-depth ratio of 4.5 using AERO-F’s DES turbulence model (based on the Spalart- 
Allmaras one-equation model (54j) and a wall-function boundary condition applied on solid surface bound¬ 
aries. The fluid domain is discretized by a mesh with 192,816 nodes and 573,840 tetrahedra (Figure[2]). The 
two-dimensional geometry is discretized in three dimensions by considering a slab of thin, but finite thickness, 
in the z-direction; the resulting grid is one element wide and is created by extruding a distance that is 1% of 
the cavity length. The viscosity is assumed to be constant, and the Reynolds number based on cavity length 
is 2.97 x 10 6 , while the free-stream Mach number is 0.6. Due to the turbulence model and three-dimensional 
domain, the number of conservation equations per node is 6, and therefore the dimension of the CFD model 
is N = 1,156,896. Roe’s scheme is employed to discretize the convective fluxes, and a linear variation of the 
solution is assumed within each control volume, which leads to a second-order space-accurate scheme on gen¬ 
eral unstructured, multi-dimensional meshes; however, we employ an extended stencil that gives fifth-order 
formal order of accuracy (with uniform mesh spacing) on inviscid, one-dimensional problems. The viscous 
flux is discretized using a centered Galerkin scheme. 



(a) Full domain 



(b) Detail around cavity 


Figure 2: Computational mesh: x — y plane cut. 


Flow simulations are performed within a time interval t G [0, T] with T = 12.5 time units. We employ 
the second-order accurate implicit three-point backward differentiation formula, which is a linear multistep 
scheme characterized by k = 2, a 0 = 1, aq = —4/3, a 2 = 1/3, (3 0 = 2/3, /3i = (3 2 = 0, for time integration; 


future work will perform numerical experiments with Runge-Kutta schemes. The OAE (2.3) arising at each 


time step is solved by a Newton-Krylov method, where GMRES is employed as the iterative linear solver 
with a restrictive additive Schwarz preconditioner (with no fill in) and the previous 50 Krylov vectors are 
employed for orthogonalization. Convergence is declared when the residual norm is reduced to a factor of 
10~ 3 of its starting value. All flow computations are performed in a non-dimensional setting. 

The initial condition Xq is provided by first computing a steady-state solution, and using that solution 
as an initial guess for an unsteady ‘transient’ simulation (which captures the initial transient before the flow 
reaches a quasi-periodic state) of 7.5 time units. The state at the end of the unsteady transient simulation is 
then used as the initial condition for the subsequent simulations. The steady-state calculation is characterized 
by the same parameters as above, except that it employs local time stepping with a maximum CFL number 
of 100, it uses the first-order implicit backward Euler time integration scheme, it assumes a linear variation 
of the solution within each control volume, it employs a Spalart-Allmaras turbulence model, and it employs 
only one Newton iteration per (pseudo) time step. 

The output of interest is the pressure at location (0.0001,-0.0508,0.0025), which is shown in the bottom 
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row of Figure |4] All errors are reported as the d 2 relative error in this quantity, i.e., 


\JyZL i‘* (i*(p)(nAt*) -p*(nAt*)) 2 

= -- 

where p : N(T/At) —>■ R is the pressure for the model of interest, p* : N(T/At*) —> K is this pressure 
response of the designated ‘truth’ model (typically the full-order model), and P* is a linear interpolation of 
the pressure response onto the grid based on the truth-model time step At*. 

All computations are performed in double-precision arithmetic on a parallel Linux clustei[^] using 48 cores 
across 6 nodes. 

7.2. Time-step verification 

Because this paper considers the time step to be an important parameter in model reduction, we first 
perform a time-step verification study to ensure we employ an appropriate ‘nominal’ time step. Figure [3] 
reports these results using a time-step refinement factor of two. A time step of At* = 0.0015 time units 
yields observed convergence rates in both the instantaneous drag force on the lower wall and instantaneous 
pressure at t = T that are close to the asymptotic rate of convergence (2.0) of three-point BDF2 scheme. 
Further, this value also leads to sub-2% errors in both quantities, which we deem to be sufficient for this set 
of experiments. 

Figure [4] shows several instantaneous snapshots of the vorticity field and corresponding pressure field 
generated by the high-fidelity CFD model. The flow within the cavity is quasi-periodic; during one cycle, 
vorticity is shed from the leading edge of the cavity, convects downstream, and impinges on the aft edge of 
the cavity. Upon impingement, an acoustic disturbance is generated which propagates upstream and scatters 
on the leading edge of the cavity, generating a new vortical disturbance to initiate the next oscillation cycle. 
The pressure fields in the bottom row of Figure [4] reveal regions of low pressure (blue contours) associated 
with vortices, as well as acoustic disturbances both within the cavity and radiating outside the cavity. This 
complex flow is governed by the interactions of several nonlinear processes, including roll-up of the shear 
layer vortices, impingement of the vortices on the aft wall resulting in sound generation, propagation of 
nonlinear acoustic waves, and interaction of these waves with the shear layer vorticity. 

7.3. Reduced-order models 

To construct both the Galerkin and LSPG ROMs, we employ the proper orthogonal decomposition (POD) 
technique; we employ a constant weighting matrix A = I for the LSPG ROM. To construct the POD basis, 
we set $ -s— & (X,v), where $ is computed via Algorithm [I] of the appendix with snapshots consisting of 
the initial-condition-centered full-order model states X = |a;*(fcA<*) — > where a;* denotes the FOM 

response computed for a time step of At* = 0.0015. Three values of the energy criterion v £ [0,1] are used 
during the experiments: v = 1 — 10 -4 (p = 204), v = 1 — 10 -5 (p = 368), and v = 1 — 10 -6 (p = 564). 
Figure [5] shows a selection of the energy component of the computed POD modes. Note that as the mode 
number increases, the modes capture finer spatial-scale behavior, which we expect to be associated with finer 
time-scale behavior; this will be verified in Section [7.5.1| 

We first repeat the time-step verification study, but we do so for the reduced-order models (again using 
the BDF2 scheme) in the time interval 0 < t < 0.55, as all Galerkin ROMs remain stable in this time interval. 
Figure [ 6 ] reports these results. First, we note that the Galerkin ROM converges an approximated rate of 2.0, 
which is what we expect given that the Galerkin ROM simply associates with a time-step-independent ODE 
( |3.2[ ). However, the LSPG ROM does not exhibit this behavior; in fact the error convergence is not even 
monotonic. This is likely due to the fact that the method does not associate with a time-step-independent 
ODE. 

We next perform simulations for both reduced-order models for all tested basis dimensions and time steps; 
Figure [7] reports the time-dependent responses. When a response stops before the end of the time interval, 




3 The cluster contains 8-core compute nodes that each contain a 2.93 GHz dual socket/quad core Nehalem X5570 processor 
with 12 GB of memory. The interconnect is a 3D torus InfiniBand. 
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(a) Drag: At* = 0.0015 yields an approximate rate of convergence of 1.94 and an estimated error in the output quantity 
(computed via Richardson extrapolation) of 1.26 x 10 —2 . The red points denote the result for At* = 0.0015. The rightmost plot 
shows the time-dependent response for the finest tested time step At = 1.875 X 10 —4 (black solid) and the converged time step 
At* = 0.0015 (red dashed). 



At 



(b) Pressure: At* = 0.0015 yields an approximate rate of convergence of 1.83 and an estimated error in the output quantity 
(computed via Richardson extrapolation) of 7.68 x 10 —4 . The red points denote the result for At* = 0.0015. The rightmost plot 
shows the time-dependent response for the finest tested time step At = 1.875 X 10“ 4 (black solid) and the converged time step 
At* = 0.0015 (red dashed). 


Figure 3: Time-step verification study. Note that the approximated convergence rates are close to the asymptotic value of 2.0 
for the BDF2 scheme. 
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(a) time = 2.10 (b) time = 2.61 (c) time = 3.12 (d) time = 3.63 



(e) time = 2.10 (f) time = 2.61 (g) time = 3.12 (h) time = 3.63 


Figure 4: Instantaneous CFD vorticity field (top) and pressure field (bottom) during one oscillation cycle. The dot on the 
forward wall of the cavity indicates the location of the pressure signal output. 


this indicates that a negative pressure was encountered, which causes AERO-F to exit the simulation. We 
interpret this phenomenon as a non-physical instability. 

First, note that the Galerkin ROMs become unstable (i.e., generate a negative pressure) for all time steps 
and all basis dimensions. This is consistent with previously reported results Gaa Bung Eg that indicate 
Galerkin projection almost always leads to inaccurate responses for compressible fluid-dynamics problems. 
In contrast, the LSPG ROM results in many stable, accurate responses for all basis dimensions. Further, 
LSPG responses exhibit a clear dependence on the time step At. Subsequent sections provide a deeper 
analysis of this dependence. 


7-4- Limiting case: comparison 

We next compare the responses of the Galerkin and LSPG ROMs for small time windows (when the 
Galerkin responses remain stable) and small time steps. Figure [8] reports e(pdiscrete opt.,PGai*) —which is the 
difference between the pressure responses generated by the LSPG ROM with different time steps and the 
Galerkin ROM with a fixed time step At = 1.875 x 10~ 4 (the smallest tested time step)—for a time window 
0 < t < 1.1. These responses support an important conclusion (see Theorem 5.3): the Galerkin and LSPG 
ROMs are equal in the limit of At —> 0 for A = l/y/aol, which is what we employ for the LSPG ROM (note 
that ao = 1 for this time integrator)]^] This has significant consequences for the LSPG ROM, as decreasing 
the time step leads to the same unstable response as Galerkin; larger time steps are needed to ensure the 
LSPG ROM is stable for the entire time interval. 

Figure [9] reports e(pdiscrete opt.>PFOM*) and £(pGai.,PFOM*) —which are the differences between the two 
ROM-generated pressure responses and the full-order model pressure response for At = 1.875 x 10 -4 —as a 
function of the time step for all three basis dimensions and three time intervals. These results highlight a 
critical observation: the LSPG ROM is more accurate for an intermediate time step. This not only supports 
the result of Corollary |6.17[ but provides an interesting insight: taking a larger time step not only leads to 
better speedups (i.e., the end of the time interval is reached in fewer time steps), but it also decreases the 
error, sometimes significantly. This is further explored in the next section. 


4 Note that in the p = 564 case, it is not clear if the difference is converging to zero. This is likely due to the fact that 
the time steps are not sufficiently small to detect convergence to zero in this case. In fact, as the basis dimension p increases, 
the basis captures finer temporal behavior (as will be shown in Figure [Tol l an<4 so the time scale of the ROM response will be 
smaller; in turn, smaller time steps At will be required to detect convergent behavior. 
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(a) mode 1 


(b) mode 21 



(c) mode 101 



(d) mode 201 


(e) mode 401 


Figure 5: Visualization of the energy component of the POD modes. 





(b) LSPG reduced-order model 


Figure 6: Time-step verification study for Galerkin and LSPG reduced-order models for p = 368 and 0 < t < 0.55. While 
the approximated convergence rate for the Galerkin reduced-order model is close to the asymptotic value of 2.0 for the BDF2 
scheme, this is not observed for the LSPG reduced-order model. Note that the error for each ROM is computed with respect its 
response for finest Richardson extrapolation, and a dashed red line indicates the snapshot-collection time step At* = 0.0015. 
See Figure [TJc—d) for the associated time-dependent responses. 


7.5. Time-step selection 

Recall from Corollary |6.17| and Remark |6.18| that decreasing the time step At has a non-obvious effect 
on the error bound for the LSPG ROM. We now assess these effects for the current problem. 


7.5.1. Spectral content of POD basis 

In our interpretation of the error bound (6.651 for the LSPG ROM applied to the backward Euler 
scheme, we noted that the time step should be ‘matched’ to the spectral content of the trial basis <I>. This 
is of practical importance, as selecting an appropriate time step for the ROM should take into account the 
relevant temporal dynamics associated with the basis. For example, a time step may be too small if the 
basis has filtered out modes with a time scale matching that of the time step. If we assume that the basis 
$ is computed via POD, then we would expect the vectors to be naturally ordered such that lower mode 
numbers are associated with lower temporal frequencies. Then, including additional modes has the effect of 
encoding information at higher frequencies. It follows that the time step should be decreased as additional 
modes are retained in construction of the ROM. 

Here we investigate the validity of this assumption by examining the spectral content of the POD basis 
vectors for the current cavity-flow problem. We compute the time histories of the generalized coordinates by 
projecting the FOM solution onto the POD basis as :r*(fcAf*) := 3> 7 (a;*(fcAf*) — a:o), k € N(8334). We then 
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(a) Galerkin, p = 204 


6 

time 

(b) LSPG, p = 204 



(c) Galerkin, p = 368 



- FOk, At = o.ool5 
A£=0.00015 
— At=0.0001875 

.Ai=0.000375 

-Af=0.00075 

-Af=0.0015 

A£=0.003 
Ai=0.006 
Ai=0.012 
Ai=0.015 



(e) Galerkin, p = 564 


6 

time 

(f) LSPG, p = 564 


Figure 7: Responses generated by Galerkin and LSPG ROMs for different basis sizes p and time steps At 
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(a) p = 204 (b) p = 368 (c) p = 564 


Figure 8: Difference £(pLSPG 5 PGai*) between the pressure responses generated by the LSPG ROM with different time steps 
and the Galerkin ROM with a fixed time step At = 1.875 x 10 —4 in 0 < t < 1.1. This demonstrates convergence of the LSPG 
ROM to Galerkin as At —> 0. 


compute power spectral densities of the generalized coordinates x+ft). Figure [l0(a)| shows sample spectra, 
normalized by the total energy in each signalj^jfor several of the POD modes. The figure shows that energy 
shifts to higher frequencies as the POD mode number increases, confirming our assumption for this example. 
This is further quantified by calculating a characteristic time-scale T 95 associated with each mode; we define 
this time scale as the inverse of the frequency below which 95 percent of the energy is captured for that 
mode. Figure [T0(b)| plots this time scale versus the mode number, showing a clear trend of decreasing time 
scale with increasing mode number. 

Thus, at least for the present application problem, we expect the optimal time step for the LSPG ROM 
to decrease as modes are added to the POD basis (this will be verified by Figure [f2|). Note that systematic 
calibration could be performed to attempt to automate selection of the ROM time step as a function of basis 
dimension. While this would be of clear practical interest, we do not pursue it here, as optimal-timestep 
computation would be complicated in practice by nonlinear interactions arising from the dynamical system, 
as well as effects from the spatial-discretization error and POD truncation error. 


7.5.2. Error bound behavior 

Having verified that higher POD mode numbers correspond to smaller wavelengths, we now numeri¬ 
cally assess quantities related to the error bound (6.65). First, Figure [Tl (a) | reports the dependence of the 
maximum relative projection error max^ ft*(<&, At) on the time step At and the basis dimension, where 


«(*,At) := 


|(7 - $$ T )(x*(fcAt) - x*((fc - 1 )At))|| 
||a;*(fcAf) - a;*((fc - 1)At)|| 


Note that ft* is closely related to ft fe from error bound ( |6.65 ), as they are equal if x a + &x P (t) = x(t) and 
the LSPG ROM computes x k P such that p k is minimized. 

These results confirm that adding basis vectors—which we know has the effect of encoding higher fre¬ 
quency content—significantly reduces the projection error for small time steps At, but has less of an effect on 
larger time steps , as re taining the first POD vectors already enables dynamics at that scale to be captured. 

Next, Figure 11(b) plots the error bound (6.65) for a value of n = 1 and with fi k = ft*. This highlights 


an important result: selecting an intermediate time step At leads to the lowest error bound, regardless of 
the basis dimension. Even though this result corresponds to the backward Euler integrator, we expect a 
similar trend to hold for the present experiment, which uses the BDF2 scheme. The next section assesses 
the performance of the LSPG ROM, including its dependence on the time step. 


5 The energy in a time series within some frequency range is obtained by integrating the power spectral density over that 
range. 
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At 

(a) 0 < t < 0.55, p = 204 



At 

(d) 0 < t < 0.55, p = 368 





(b) 0 < t < 1.1, p = 204 (c) 0 < t < 1.54, p = 204 




At 


At 


(e) 0 < t < 1.1, p = 368 



(f) 0 < t < 1.65, p = 368 



At 


At 


At 


(g) 0 < t < 0.55, p = 564 (h) 0 < t < 1.1, p = 564 (i) 0 < t < 1.65, p = 564 

Figure 9: Galerkin errors e(pGal., PFOM* ) and LSPG errors e(pLSPGi PFOM* ) over different time intervals, time steps, and basis 
dimensions. For reference, a dashed red line indicates the snapshot-collection time step At* = 0.0015. 
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(a) Power spectral densities of several generalized coordi- (b) Dependence of the POD mode characteristic time scale 
nate time series for the cavity-flow problem. on mode number for the cavity-flow problem. 


Figure 10: Spectral content of the POD basis. 


7.6. LSPG ROM performance 

We now compare the accuracy and walltime performance of the LSPG ROM as the dimension of the 
basis, time step, and time interval change. The most salient result from Figure [12] is that choosing an 
intermediate time step leads to both better accuracy and faster simulation times. This shows that our 
theoretical analysis of the error bound performed in Section |7.5.2| leads to an actual observed performance 
improvement. For example, consider the p = 564 case over the time interval 0 < t < 2.5. In this case, 
a time step of At = 1.875 x 10~ 4 leads to a relative error of 0.0140 and a simulation time of 289 hours; 
increasing this value to At = 1.5 x 10 -3 reduces the relative error to 9.46 x 10 -4 and the simulation time 
to 35.8 hours, which constitutes roughly an order of magnitude improvement in both quantities. Again, this 
supports the theoretical results of Corollary |6 .1 7| and highlights the critical importance of the time step for 
LSPG reduced-order models. 

In addition, Figure [12] shows that as the basis dimension increases, the optimal time step decreases; this 
was anticipated from the spectral analysis performed in Section [7.5.1| In addition, adding POD basis vectors 
does not improve accuracy for large time steps. We interpret this effect as follows: for larger time steps, 
the first few POD modes accurately capture ‘coarse’ phenomena on the scale of the time step. Therefore, 
accuracy improvement is not achieved by adding modes that encode dynamics that evolve on a time scale 
finer than the time step itself. 

Further, Figure [12(g)] highlights that as the basis dimension increases, the error generally decreases, which 
is an artifact of the monotonic decrease in the FOM OAE residual achieved by the LSPG ROM (Remark 
4.1). Finally, the figure shows that as the time interval grows, the optimal time step generally increases. 


7.7. GNAT: ROM with complexity reduction 

In this section, we perform a similar study, but equip the LSPG ROM with complexity reduction in order 
to achieve computational savings. In particular, we employ the GNAT method mmm, which solves 
Eq. (4.1) with A = (P«I> r ) + P, where <J> r is a basis for the residual and P consisting of selected rows of the 


identity matrix. 

The problem is identical to that described in Section |7.1| except that we take T = 5.5 time units and 
employ a second-order space-accurate dissipation scheme wherein a linear variation of the solution is assumed 
within each control volume]^] For this simulation, the full-order model consumes 5.0 hours on 48 cores across 
six compute nodes. 

To construct the trial basis $ and basis for the residual 3> r for the GNAT models, we again employ POD. 
In particular, we set $ •<— $ (A, zz), where $ is computed via Algorithm [I] with snapshots consisting of the 


3 This is done to ensure the sample mesh requires two layers of neighboring nodes for each sample node. 
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(a) Dependence of the maximum relative projection error (b) L SPG error bound for the backward Euler scheme 
maxfc fl k on the time step At and basis dimension p (|6.65| using experimental data for fi k and setting k, = 1 

and p k = 


Figure 11: Assessment of quantities appearing in error bound ( |6.65| >. This analysis suggests that an intermediate time step At 
can reduce errors for the LSPG ROM. 


centered full-order model states X = {:r*(fcAt*) — ■ An energy criterion of v = 1 — 10 5 (p = 179) is 

used during the experiments. For the residual, we employ «— <I> (X r . v r ) via Algorithm [l] with snapshots 
X r = {r n (x o + 3>w n< ' k ' > ), k £ N(K(n)), n £ N(2228)} and w n ^ corresponding to the LSPG ROM solution 
at Gauss-Newton iteration k within time step n using a time step of At = 6 x 10 -3 . Here, K(n) denotes the 
number of Newton iterations required for convergence of at time instance n. An energy criterion of v r = 1.0 
is employed. In addition, the GNAT model sets the Jacobian basis equal to residual basis 3>,/ = 3> r and 
employs n s = 743 sample nodes that define P, which leads to 4458 rows in P as there are six conservation 
equations per node due to the turbulence model (see Ref. [2T] for definitions). 

The GNAT implementation in AERO-F is characterized by the sample-mesh concept [211; . Figure [l3| 
depicts the sample mesh for this problem, which was constructed using n c = 2228 working columns | 21 [ 
Algorithm 3], and includes two layers of nodes around the sample nodes (to enable the residual to be 
computed at the sample nodes). It is characterized by 7,974 total nodes (4.1% of the original mesh) and 
17,070 total volumes (3.0% of the original mesh). Due to the small footprint of the sample mesh, the GNAT 
simulations are run using only 2 cores on a single compute node. 

Figure [14] reports the results obtained with the GNAT ROM using different time steps. Critically, note 
that the GNAT ROM also exhibits a ‘dip’ in the optimal time step, with a time step of 6.0 x ICR 3 yielding 
the lowest error. In fact, increasing the time step from 1.5 x 10 -3 to 6.0 x 10 -3 decreases the error from 
3.32% to 2.25% and also significantly increases the computational savings relative to the full-order model 
(as measured in core-hours) from 14.9 to 55.7. This highlights that the analysis is also relevant to ROMs 
equipped with complexity reduction. 

7. 8. Summary of experimental results 

We now briefly summarize the main experimental results: 

• Galerkin ROMs are unstable for long time intervals (Figure [7|. 

• LSPG ROMs are only unstable for small time steps (Figure [7]). 

• Galerkin and LSPG ROMs are equivalent as At —>• 0 (Figure [ 8 |. 

• LSPG ROMs are more accurate than Galerkin ROMs over small time windows where Galerkin is stable 
(Figure [9|. 

• LSPG ROMs are most accurate for an intermediate time step (Figure [9]). 
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(g) Dependence of optimal time step (in terms of minimizing 
error) on time interval and basis dimension 


Figure 12: Dependence of error and simulation time for the LSPG reduced-order model on the time step At, basis dimension, 
and time interval 
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(a) Full domain 



(b) Zoom on cavity 


Figure 13: Sample mesh (red) embedded within original mesh. The sample mesh is used to to enable low-computational-footprint 
simulations with the GNAT reduced-order model. 



time 


(a) responses 



(b) errors 



Figure 14: Responses, errors £(pgnat?Pfom* )? and computational savings (as measured in core-hours) produced by the GNAT 
reduced-order model for different time steps At. 
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• Adding POD modes has the effect of including higher-frequency response components (Figure [To|. 

• The theoretical error bound for the LSPG ROM exhibits the same time step ‘dip’ as the experimentally 
observed error (Figure [TTj) . 

• The optimal time step for the LSPG ROM decreases as modes are added to the POD basis (Figure 

121 . 

• Adding modes to the POD basis has little effect on LSPG ROM accuracy for large time steps (Figure 
U). 

• The optimal time step for the LSPG ROM tends to increase as the time interval increases (Figure 

BH >- 

• The GNAT ROM, which is discrete optimal and is equipped with complexity reduction, also produces 
minimal error for an intermediate time step (Figure |l4|). 


8. Conclusions 

This work has performed a comparative theoretical and experimental analysis of Galerkin and LSPG 
reduced-order models for linear multistep schemes and Runge-Kutta schemes. We have demonstrated a 
number of new findings that have important practical implications, including conditions under which the 
LSPG ROM has a time-continuous representation, conditions under which the two techniques are equivalent, 
and time-discrete error bounds for the two approaches. 

Perhaps most surprisingly, we demonstrated that decreasing the time step does not necessarily decrease 
the error for the LSPG ROM. This phenomenon arose in both the theoretical analysis and in numerical 
experiments. In particular, our results suggest that the time step should be ‘matched’ to the spectral content 
of the reduced basis. In the experiments, we showed that increasing the time step to an intermediate value 
decreased both the error and the simulation time by an order of magnitude in certain cases. Alternatively, 
decreasing the time step cause the LSPG ROM to become unstable for longer time intervals. This highlights 
the critical importance of time-step selection for LSPG ROMs. 
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Appendix 

Algorithm [T] reports the algorithm for computing a POD basis using normalized snapshots. 
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Input: Set of snapshots X = [wijfjtj C energy criterion v G [0,1] 

Output: $ (X, v) 

1 : Compute thin singular value decomposition W = UT,V T , where W = [mi/||iui|| • • • «Jn„/||^n„||] • 
2 : Choose dimension of truncated basis p = n e (v), where 

n e (v) = arg min i 
' 8 ieV(v) 

n n w 

V( v ) := {n G N(n w ) | ^ erf/^ of > v}, 

i= 1 i=1 


and S = diag (u,;). 

3: $ (X , is) = [w 1 • • • m p ] , where U = [u 1 • • • u 71 ™]. 

Algorithm 1: Proper-orthogonal-decomposition basis computation (normalized snapshots) 
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